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This  work  examines  the  feasibility  of  and  requirements  for  using  a Shack- 
Hartmann  wavefront  sensor-based  system  for  the  tomographic  reconstruction  of  fluid 
density  in  an  aerodynamic  flow  measurement.  The  sampling  and  noise  requirements  are 
examined  at  each  stage  of  the  measurement  process,  including  the  Shack-Hartmann 
sensor,  phase  reconstruction,  and  tomographic  reconstruction.  Based  on  this  analysis,  the 
key  design  relationships  for  the  Shack-Hartmann  sensor  are  developed.  A two- 
dimensional  Shack-Hartmann  sensor  was  designed  and  built  using  a commercial  off-the- 
shelf  lens  array  and  an  8-bit  CCD  camera.  The  characterization  and  calibration  test 
results  are  presented.  The  system  capabilities  are  demonstrated  by  taking  a single 
projection  measurement  of  a free  air  jet,  and  a simulated  tomographic  reconstruction  is 
performed  by  assuming  axial  symmetry.  The  results  indicate  that  the  Shack-Hartmann 
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sensor  is  a viable  tool  for  this  type  of  measurement  and  that  it  may  be  especially  well 
suited  to  the  transonic  and  subsonic  measurement  regime. 
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CHAPTER  1 
INTRODUCTION 

Problem  Statement 

In  the  study  of  aerodynamic  flow  patterns,  most  fluids  are  transparent  media  so  that 
their  motion  remains  invisible  to  the  human  eye.  Significant  insight  is  gained  by  making 
this  flow  pattern  visible.  These  flow  visualization  techniques  have  played  an  important 
role  in  the  study  of  fluid  dynamics  (Rowley  1969;  Matulka  & Collins  1971).  In  addition 
to  providing  a qualitative  visualization  of  the  flow  phenomena,  many  of  these  techniques 
are  capable  of  providing  quantitative  data  (Rowley  1969;  Zien  et  al.  1974).  The  most 
valuable  flow  visualization  methods  are  those  that  are  able  to  supply  quantitative  data 
over  the  entire  region  of  interest  without  interfering  with  the  flow.  Optical  methods, 
which  are  sensitive  to  changes  in  the  refractive  index  of  the  flow  field,  are  especially  well 
suited  for  this  type  of  non-intrusive  measurement.  The  approach  taken  in  this  work  will 
be  to  examine  the  feasibility,  requirements,  and  anticipated  performance  of  a Shack- 
Hartmann  wavefront  sensor  based  system  for  performing  the  type  of  quantitative  flow 
measurement  needed  for  tomographic  reconstruction. 

In  order  to  reconstruct  a three-dimensional  density  profile,  the  principles  of 
tomographic  reconstruction  can  be  used  to  find  the  two-dimensional  profiles  for  multiple 
slices  of  the  test  region,  which  are  stacked  to  form  the  three-dimensional  reconstruction. 
Projection  data  must  be  available  over  a full  180  degree  field  of  view,  surrounding  each 
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'slice',  to  implement  this  approach  (Kak  & Slaney  1988).  In  an  optical  tomographic 
reconstruction,  the  variations  in  refractive  index  of  the  test  region  are  reconstructed.  The 
variations  in  refractive  index  are  manifested  in  the  form  of  optical  phase  shifts  that  can  be 
related  to  the  fluid  density.  Figure  1-1  shows  a conceptual  block  diagram  of  the 
measurement  of  projection  data  for  a single  projection  angle.  The  source  provides  a 
plane  wave  reference  beam,  which  passes  through  the  test  region.  The  wavefront  sensor 
measures  the  spatial  phase  distribution  of  the  wavefront  after  it  has  passed  through  the 


Figure  1-1:  Optical  projection  data 


test  region.  In  order  to  apply  tomography  to  reconstruct  a two-dimensional  slice  in  the 
y-z  plane,  multiple  measurements  must  be  made  for  varying  projection  angles,  y, 
covering  the  full  field  of  view  (y  = 0 tol  80  degrees). 

Research  Objectives 

In  order  for  a wavefront  sensor  to  be  used  in  a tomographic  flow  visualization 
system,  the  tomographic  reconstruction  requirements  and  constraints,  as  well  as  the  test 
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object  characteristics,  must  be  tied  to  the  wavefront  sensor  design.  Figure  1-2  shows  how 
the  various  requirements  and  characteristics  are  linked  to  the  wavefront  sensor  design. 


Tomographic 

Reconstruction  Requirements 


> Reconstruction  Resolution 

> Number  of  Projection  Angles 

> Projection  Data  Sampling 

> System  Constraints 


Test  Object 
Characteristics 


> Spatial  Fretpiency 

> Magnitude  of  Phase  Shift 
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WaveFront  Sensor 
Design 

> Lenslet  Aperture  Size 

> Lenslet  Shape 

> Lenslet  Spacing 
>Lens  Array  Size 

> Lenslet  Focal  Length 

> Focal  Spot  Centroiding 

> CCD  Geometry 
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Figure  1-2:  Wavefront  sensor  design  trade-offs/interdependencies 


The  goal  of  the  present  effort  is  to  examine  and  develop  the  theory  required  to  build 
a two-dimensional  Hartmann  wavefront  sensor  system  for  application  in  general  flow 
visualization  and  to  experimentally  verify  the  measurement  capabilities  of  the  system. 
The  relationships  and  the  interactions  between  expected  wavefront  characteristics,  lens 
array  requirements,  CCD  array  requirements,  spot  deviation  measurement  accuracy  and 
noise  characteristics,  and  tomographic  reconstruction  all  need  careful  examination  if  they 
are  to  be  interconnected  to  form  an  effective  measurement  tool. 

The  primary  analytical  result  of  this  work  will  be  to  outline  the  design  procedure  for 
a Hartmann  sensor  based  system  to  be  used  for  projection  measurement  in  tomographic 
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reconstruction.  Given  an  available  camera,  or  CCD  array,  and  reconstruction  resolution 
requirement,  the  process  of  designing  all  of  the  optics,  including  the  lens  array  will  be 
carried  out.  In  order  to  tie  all  of  the  analytical  results  together  and  to  the  overall  research 
goal,  a simple  design  example,  which  is  based  on  a supersonic  conical  flow  problem,  will 
be  carried  through  each  step.  Each  chapter  will  conclude  with  a summary  including  the 
application  of  the  theory  to  the  design  example.  The  analytical  tasks  can  be  outlined  as 
follows. 

Projection  Data  Format 

A determination  of  the  required  number  of  projection  data  points  for  a given 
tomographic  resolution  requirement  is  one  of  the  first  steps  in  the  process.  This  will  be 
dependent  upon  the  spatial  frequency  characteristics  of  the  flow  field  and  the  number  of 
projection  angles  available. 

Projection  Data  Noise  Requirements 

The  overall  accuracy  of  the  tomographic  reconstruction  will  be  partially  dependent 
upon  the  accuracy  of  the  projection  data.  The  ability  to  place  a quantitative  requirement 
on  the  quality  of  the  projection  data  will  therefore  be  crucial  to  the  design  of  the 
wavefront  sensor.  To  accomplish  this,  a statistical  analysis  must  be  performed  on  the 
tomographic  reconstruction  process.  This  analysis  should  be  tractable  for  the 
deterministic  approaches  such  as  Fourier  transform  and  filtered  back-projection 
reconstruction.  For  the  iterative  approaches,  there  is  a broad  base  of  knowledge  in  the 
literature  concerning  this  requirement  which  has  been  obtained  by  Monte-Carlo 
simulation  analysis. 
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Projection  Data  Generation/Wavefront  Reconstruction 

Given  the  requirements  established  above,  an  approach  to  the  generation  of  the  actual 
projection  data  must  be  established.  The  approaches  could  span  a broad  range  of 
complexities.  The  simplest  approaches  could  be  a direct  translation  of  spot  deviation  to 
phase  slope,  which  could  be  integrated  discretely  across  the  aperture  to  obtain  phase. 
More  complex  approaches  would  be  least  squares  fitting  and  modal  expansions.  There  is 
a significant  amount  of  information  from  the  literature  regarding  this  problem  as  related 
to  atmospheric  distortion.  A variation  of  one  or  more  of  these  methods  will  be  used  in 
this  work. 

Wavefront  Sensor  Design 

The  design  of  the  wavefront  sensor  must  take  into  account  all  of  the  requirements 
from  the  above  analysis,  as  well  as  the  characteristics  of  the  measurement  field,  and  the 
practical  availability  of  critical  components.  The  lens  array,  used  in  the  hartmann  sensor, 
will  be  restricted  by  the  state  of  the  art  in  microlens  array  manufacturing,  and  this  will 
result  in  some  system  design  limitations.  The  measurement  of  spot  deviation  on  the  CCD 
camera  is  the  basis  upon  which  the  system  operates.  An  analysis  of  the  spot  deviation 
dynamic  range  and  minimum  sensitivity  is  required  to  determine  the  number  and 
geometry  of  CCD  pixels  for  each  lens  array  sub-aperture.  A noise  analysis  of  this 
measurement  is  also  required  to  establish  the  measurement  accuracy. 

Design  Example 

The  design  example  to  be  used  for  this  work  will  be  that  of  a supersonic  conical 
projectile.  A numerical  model  for  a mach  1.5  conical  flow  will  be  used  to  determine  the 
characteristics  of  the  flow  field.  The  test  region  of  interest  will  consist  of  a 50mm  cube 
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containing  the  tip  of  the  conical  body,  and  the  minimum  resolvable  feature  size  will  be 
1mm.  Refer  to  Figure  1-3. 

The  design  example,  which  will  be  carried  through  each  stage  of  analysis  and  design 
is  based  on  a fluid  density  measurement  of  a mach  1 .5  supersonic  flow  over  a right 


circular  conical  body  at  zero  angle  of  attack.  This  measurement  will  be  used  as  the  basis 
for  the  design  of  the  final  system.  Due  to  the  formulation  of  the  aerodynamic  model, 
which  is  discussed  in  the  next  chapter,  the  cone  angle  will  be  chosen  such  that  a 60° 
oblique  shock  is  produced.  The  design  example  parameters  are  summarized  in  Table  1-1 . 

Experimental  Verification 

To  verify  the  analytical  results  of  this  work,  a prototype  system  capable  of 
measuring  two  dimensional  projection  data  will  be  built  and  tested.  For  the  purposes  of 
this  work,  a standard  monochromatic  CCD  camera  will  be  used.  Since  multiple  pixels 
are  required  for  each  Hartmann  sub-aperture,  the  camera  must  have  high  pixel  density 
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Table  1-1:  Design  example  parameters 


Flow  Velocity 

Mach  1.5 

Shock  Angle 

o^ 

o 

o 

Fluid 

Air  at  standard  conditions 

Test  Region  Size 

50mm  cube 

Reconstruction  resolution 

<lmm 

Optical  Wavelength 

532nm 

(512  X 512  or  higher)  in  addition  to  readily  available  access  to  the  digital  data.  Due  to  the 
anticipated  high  cost  of  a custom  lens  array,  it  is  desired  that  a commercially  available 
lens  array  be  used.  As  a result,  the  available  lenslet  aperture  sizes  will  be  a primary 
factor  in  the  determination  of  sensor  resolution. 

Static  Testing 

To  verify  the  accuracy,  sensitivity,  dynamic  range,  and  repeatability  of  the  wavefront 
reconstruction  from  the  hartmann  sensor  system,  static  tests  will  be  run  allowing  the 
measured  data  to  be  compared  to  known  data.  A simple  tilting  of  the  wavefront  can  be 
accomplished  by  a mechanically  adjusted  fold  mirror.  By  sampling  a portion  of  the 
beam,  prior  to  the  wavefront  sensor,  an  external  wavefront  tilt  measurement  can  be  done 
as  a reference  for  the  measurement. 

Dynamic  Testing 

Due  to  limitations  in  available  equipment,  dynamic  testing  will  be  very  limited. 
Because  of  the  lack  of  an  externally  triggered  frame  grabber,  any  testing  will  have  to  be 
done  using  a steady  state  flow.  Consequently,  a free  air  jet  nozzle  will  be  used  as  a 
simple  demonstration  test.  Note  that  the  exact  density  structure  of  the  nozzle  flow  will  be 
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impossible  to  establish  and  quantify  using  the  simple  lab  setup  available.  Also,  the 
supersonic  nozzle  flow  will  contain  numerous  shocks,  for  which  the  Shack-Hartmann 
sensor  is  not  well  suited.  However,  since  the  aerodynamic  community  is  very  familiar 
with  the  jet  nozzle  flow  problem,  it  is  anticipated  that  the  inclusion  of  this  demonstration 
test  will  provide  additional  insight. 

The  next  chapter  will  provide  some  background  on  the  Shack-Hartmann  sensor  as 
well  as  discuss  the  relevant  literature  for  the  general  flow  measurement  application.  This 
will  be  followed  by  a discussion  of  tomographic  reconstruction  in  chapter  three.  Chapter 
four  outlines  the  phase  reconstruction  process,  which  is  used  to  determine  the  wavefront 
phase  from  the  Shack-Hartmann  measurements,  and  chapter  five  details  the  spot 
deflection  measurement  operation  of  the  wavefront  sensor.  With  the  theory  developed  in 
these  chapters,  chapter  six  will  outline  the  wavefront  sensor  design  approach.  Chapter 
seven  will  discuss  the  measurement  results,  and  chapter  eight  will  offer  conclusions  and 


recommendations  for  future  research. 


CHAPTER  2 
BACKGROUND 


Gladstone-Pale  Relation 

In  many  applications,  we  are  interested  in  the  density  of  a gas  contained  in  the  region 
surrounding  the  aerodynamic  event  of  interest.  Under  the  assumption  that  the  material  is 
a gas  and  that  the  wavelength  of  light  being  used,  X , is  far  from  the  resonant  wavelength 
of  the  gas,  X- , the  harmonic  electron  oscillator  model  (Lorentz  Model)  can  be  used  to 
relate  the  polarization  vector  to  the  electric  field  vector.  This  will  generate  the  Claussius- 
Mosotti  relation  (Merzkirch  1987) 


.V. 


-1 


n^+2  3Km„M^v]-v^ 


[2.1] 


which  is  expressed  in  terms  of  refractive  index  n , fluid  density  p , Loschmidts  constant 
92 , electron  charge  e , electron  mass  , molar  weight  M,  molecular  oscillator  strengths 
f. , molecular  resonant  frequencies  v- , and  light  frequency  v . With  the  refractive  index 
of  most  gases  being  close  to  one  the  - 1 term  can  be  approximated  as  2(n  - 1)  and  the 
+2  term  can  be  approximated  as  3,  resulting  in  tbe  Gladstone-Dale  relation. 


n-1 


f- 


l/rm^  M i v] 


or  n-\-  Kp  . 


[2.2] 


9 


10 


K is  the  Gladstone-Dale  constant  and  is  dependent  on  the  characteristics  of  the  gas  as 
well  as  the  wavelength  of  light  used.  Values  of  K for  various  gases  can  be  found  in  many 
textbooks  and  handbooks  on  physical  chemistry.  For  air  at  a temperature  of  288K,  the 


Gladstone-Dale  constant  ranges  from  0.2239  at  a wavelength  of  0.9125  jjm  to 
0.2330  at  a wavelength  of  0.3562  jum  showing  that  it  is  only  weakly  dispersive  for 
wavelengths  far  from  resonance  (Merzkirch  1987). 

The  dynamic  range  and  spatial  frequency  content  of  the  density  field  will  be 
dependent  upon  the  type  of  flow  measurement  being  performed.  Supersonic  flow 
measurements  and  others  containing  shock  waves  will  constitute  an  upper  limit  on  spatial 
frequency.  Here  the  density  change  across  the  shock  front  is  approximately  a step 
function  since  the  change  occurs  over  a few  microns.  The  magnitude  of  the  density 
change,  for  non-isentropic  and  adiabatic  conditions  across  an  oblique  incidence  shock, 
can  be  closely  approximated  by  the  following  (Clancy  1975): 

p,  _y  + l M^sin^(g)  3] 

^ l + ^^M'sin'(f) 

2 

In  [2.3]  Y is  the  adiabatic  index  (-1.4  for  air),  M is  the  Mach  number,  s is  the  shock 

incidence  angle,  and  ^ is  the  ratio  of  density  behind  the  shock  to  that  in  front  of  it.  This 

relationship  shows  that  the  ratio  of  densities  across  the  shock  approaches  a limit  of  6 for 
normal  incidence  and  high  Mach  number.  Equation  [2.3]  can  be  used  to  relate  this 
density  ratio  to  a ratio  of  refractive  index. 

Supersonic  Conical  Flow  Model 

A numerical  solution  to  the  supersonic  conical  flow  problem  was  developed  using 
the  Taylor  and  Maccoll  solution,  as  outlined  by  Anderson  (1990),  using  Matlab.  In  this 
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solution,  a shock  angle  is  first  assumed  and  the  flow  conditions  across  the  oblique  shock 
are  determined  assuming  steady  adiabatic  flow.  Figure  2-1  Shows  the  geometry  used  for 
this  solution. 


The  differential  equations  governing  the  flow  between  the  shock  and  the  conical 
body  for  axisymmetric,  adiabatic,  steady  flow  are  then  solved  numerically  starting  from 
the  initial  condition  at  the  shock  boundary  and  progressing  toward  the  body.  The 
assumption  of  an  infinite  conical  body  results  in  a solution  that  is  only  dependent  upon 
the  angle,  0,  to  the  field  point.  At  the  conical  body,  the  normal  flow  in  the  6 direction 
will  go  to  zero  at  which  point  the  numerical  solution  is  halted. 

At  the  shock  boundary  the  free  stream  mach  number,  , can  be  separated  into 
tangential  and  normal  components,  Mti  and  Mni,  where 


V 


00 


> 


Figure  2-1:  Supersonic  conical  flow 


M,x  = cos  and 
= M„sin6>,  . 


[2.4] 
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Assuming  steady  flow  and  no  body  forces,  the  tangential  component  of  the  flow 
velocity  is  continuous  across  the  shock  boundary  (Anderson  1990).  For  a steady, 
adiabatic  flow  with  no  body  forces,  the  normal  velocity  component  just  inside  the  shock 
boundary  can  be  found  as  follows  (Anderson  1990). 


^2  = 


<1  + 

1 

1 

— 1 
1 

1 1 

A 

1 I 

Ml,  - 1 

[2.5] 


Again,  y is  the  adiabatic  index  of  the  fluid  and  is  approximately  1 .4  for  air  at 
standard  conditions.  Now,  using  a non-dimensional  velocity  defined  as  follows. 


V = 


'(y-l)M' 


+ 1 


the  governing  differential  equations  become. 


[2.6] 
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[2.7] 


= 0. 


The  above  differential  equations  were  solved  numerically  using  a 4**’  order  Runge- 
Kutta  method  beginning  from  the  initial  conditions  just  behind  the  shock  and  moving  in 
uniform  increments  of  the  angle,  0,  toward  the  body  surface  (Jain  et  ah,  1985).  The  result 
of  the  above  solution  is  a known  fluid  velocity  for  any  point  between  the  shock  and  the 
body  surface,  which  is  dependent  only  upon  the  angle,  0.  The  fluid  density  can  then  be 
found  for  isentropic  conditions  using  the  following  relationship  (Anderson  1990): 
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1 + ^^M" 
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[2.8] 
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The  results  of  the  numerical  solution  were  verified  by  comparison  of  several  test 
cases  to  the  results  published  in  tables  from  Sims  (1964).  The  Gladstone-Dale  relation 
can  now  be  used  to  relate  the  fluid  density  to  the  index  of  refraction.  Figure  2-2  shows  a 
plot  of  the  index  profile  as  a function  of  the  radius  from  the  axis  of  the  conical  body. 


Figure  2-2:  Refractive  index  profile 


With  the  refractive  index  profile  known  for  all  space,  the  phase  shift  of  a plane  optical 
wave  passing  through  that  space  can  be  found.  Figure  2-3  shows  the  geometry  used  for 
the  solution  of  this  problem  for  a cross  sectional  slice  in  the  x-y  plane,  located  at  z=Zq. 
Outside  the  shock  boundary  the  index  of  refraction  will  be  constant  while  the  region 
between  the  shock  and  the  body  will  have  an  index  that  varies  with  the  angle,  0.  For  a 
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plane  wave  incident  along  the  y direction  and  normal  to  the  z-axis,  there  will  be  a phase 
shift  relative  to  a portion  of  the  beam  outside  the  shock  region,  which  is  a function  of  x 
and  can  be  written  as  follows: 

ry  y„ 

^<l>  = ^ \[n{x,y,z^)-n^]dy  . [2.9] 

->■» 


Figure  2-3:  Geometry  and  coordinate  system 


In  equation  [2.9]  X.  is  the  wavelength  and  the  integration  limits  +yo  can  be  written  as  a 
function  of  x,  as  follows: 


[2.10] 


Note  that  the  phase  for  |x|  < is  undefined  due  to  the  obscuration  of  the  conical  body 
and  should  be  set  to  zero.  Since  the  index  profile  is  dependent  only  upon  the  angle  0,  the 
above  calculations  can  be  done  for  a normalized  case  of  rs  =1,  and  scaled  for  any  other 


15 


value  of  Ts.  This  integrated  phase  shift  constitutes  one  projection  angle  of  the  projection 
data  required  for  tomographic  reconstruction  of  a slice  of  the  index  profile  at  z=Zq.  In 
this  axisymmetric  case,  no  other  projection  angles  are  required  since  there  is  no 
dependence  on  rotation  angle.  Figure  2-4  shows  an  example  plot  of  the  integrated  phase 
shift  (projection  data)  for  a mach  1 .5  flow  with  a shock  angle  of  60  degrees. 


Figure  2-4:  Integrated  phase  (mach=1.5,^^  =6Q°,r^  = \ cm,  A = 532  nm) 

Note  that  the  phase  slope  approaches  infinity  as  the  shock  is  approached  as  expected. 
In  order  to  build  up  the  entire  three-dimensional  index  profile  of  the  region  between  the 
shock  and  the  body,  multiple  slices  for  varying  values  of  Zo  are  required  so  that  the  above 
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calculations  are  required  for  all  points  in  the  x-z  plane.  Figure  2-5  shows  this  phase  data 
in  an  interferogram  format  for  the  same  test  conditions  as  above.  In  this  type  of  plot  the 
intensity  that  would  be  recorded  by  the  interference  of  a beam  passing  through  the  test 
region  and  a plane  reference  wave  is  shown  in  the  interferogram.  The  results  of  the 
above  analysis  will  be  applied  to  the  design  example  at  the  end  of  the  chapter  in  order  to 
determine  the  maximum  allowable  phase  slope  for  the  projection  measurements. 

Non-Intrusive  Flow  Measurement  Techniques 
Shadowgraph  and  Schlieren  Photography 

Shadowgraph  and  Schlieren  methods  were  among  the  earliest  flow  visualization 
methods  used.  Qualitative  visualization  is  obtained  by  these  methods  which  are  sensitive 
to  changes  in  the  first  and  second  derivatives  of  the  density.  Although  Weinstein 
developed  a modified  focused  schlieren  system,  which  was  capable  of  making 
quantitative  measurements  in  a focussed  plane  within  the  flow  field,  these  methods  are 
generally  not  suitable  for  obtaining  quantitative  data  (Weinstein  1991). 

Interferometry 

Optical  Interferometry  has  long  been  used  to  reconstruct  the  refractive  index  field  of 
a transparent  medium  by  measuring  the  phase  shift  of  light  rays  passing  through  the 
medium  (Bockasten  1961;  Barr  1962;  Maldonado  et  al.,  1965;  Maldonado  & Olsen  1966; 
Olsen  et  al.  1968;  Bruning  et  al.  1974).  By  recording  the  interference  pattern  obtained  by 
the  superposition  of  a test  wave  passing  through  the  medium  of  interest  and  a separate 
reference  wave,  the  optical  path  differences  can  be  obtained.  Figure  2-6  shows  an 
example  of  how  this  measurement  can  be  made  with  a Mach-Zehnder  interferometer. 
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Figure  2-5:  Phase  interferogram  2-D  projection  data 
(Mach=l  .5, = 60° , /I  = 532  nm ) 

The  measurement  is  performed  by  recording  the  intensity  of  the  interference  pattern 
in  the  output  recording  plane  which  will  produce  a pattern  similar  to  Figure  2-5.  Given 
the  recorded  intensity  values,  the  straightforward  relationship  between  the  phase  of  the 
test  beam  and  the  interferogram  intensity  can  be  used  to  solve  for  the  test  beam  phase. 


18 


The  optical  path  difference  represents  the  integrated  phase  shift  incurred  in  passing 
through  the  test  medium  due  to  the  inhomogeneous  refractive  index.  For  the 
axisymmetric  case,  the  Abel  transform  can  be  used  to  invert  the  integrated  phase  and 
reconstruct  the  refractive  index  profile.  In  the  general  case,  where  the  field  is  not  axially 
symmetric,  multidirectional  interferometric  data  is  needed  and  the  principles  of 
tomography  can  be  applied  to  reconstruct  the  three-dimensional  refractive  index  field. 


Mirror 


Fig  2-6:  Interferometry  using  a Mach-Zehnder  interferometer 


Holographic  Interferometry 

Soon  after  the  introduction  of  off-axis  holography  by  Leith  and  Upatnieks  (1962) 
this  new  technique  was  looked  upon  as  a potentially  useful  tool  in  the  field  of 
interferometry.  As  demonstrated  in  the  work  by  Leith  and  Upatnieks,  if  a test  (object) 
beam  is  superimposed  with  an  off  axis  reference  beam  and  the  interference  pattern  is 
recorded  on  film,  then  the  test  beam  can  be  reconstructed  later  using  the  developed  film, 
(hologram),  and  only  the  off-axis  reference  beam.  Figure  2-7  shows  the  process  of 
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Figure  2-7:  Holography,  (a)  recording  and  (b)  reconstruction 


recording  the  hologram  and  using  it  to  reconstruct  the  recorded  object  wave. 

In  interferometry  this  would  allow  a hologram  to  be  placed  in  the  test  leg  of  an 
interferometer  instead  of  the  actual  test  field  as  first  pointed  out  by  Horman  (1965). 

Since  the  holographic  technique  was  capable  of  essentially  freezing  a transient  event  in 
time  by  capturing  it  on  film,  this  represented  a significant  extension  of  the  field  of 
interferometry  into  the  study  of  aerodynamics.  At  around  this  same  time  Stetson  and 
Powell  (1965)  demonstrated  the  use  of  holographic  interferometry  in  real  time  vibration 
analysis. 

Theory  of  operation 

Holographic  interferometry  has  been  used  extensively  in  flow  visualization  (Horman 
1965;  Heflinger  et  al.l966;  Rowley  1969;  Matulka  & Collins  1971;  Zien  et  al.  1974). 
Here,  a hologram  is  recorded  of  the  flow  field,  which  allows  the  interference  pattern  to  be 
reconstructed  later  using  a reference  wave  only.  The  off-line  reconstruction  simplifies 
the  recording  of  data  and  allows  it  to  be  taken  over  a limited  field  of  view  (+  5 to  10 
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degrees).  Multiple  holograms  generally  taken  over  a 180-degree  field  of  view  are 
required  to  reconstruct  the  flow  field. 

Holographic  interferometry  for  fluid  flow  measurement 

The  first  direct  demonstration  of  the  use  of  holographic  interferometry  in  the  study 
of  high-speed  aerodynamic  events  was  published  by  Heflinger  et  al.  in  February  1966.  In 
this  work  the  authors  have  done  away  with  the  conventional  interferometer  unlike  the 
idea  put  forth  by  Horman  (1965).  Here  they  were  interested  in  reconstructing  the 
interference  pattern  between  the  reference  and  test  waves  instead  of  reconstructing  the 
test  wave  itself  To  do  this,  a double  exposure  was  made  in  two  separate  steps  resulting 
in  an  interferogram.  The  first  exposure  is  made  just  prior  to  the  flow  event  taking  place 
representing  the  ‘comparison  scene’  and  the  second  one  completes  the  exposure  at  the 
instant  desired  during  the  flow  event  representing  the  ‘test  scene’.  The  process  is  shown 
in  Figure  2-8.  The  reference  beam  for  both  exposures  is  identical  except  in  absolute 
phase.  After  the  interferogram  is  developed  the  reference  beam  is  used  with  it  to 
reconstruct  both  ‘scenes’  by  superposition  which  forms  the  interference  ‘fringe’  pattern 
representing  the  difference  between  the  test  scene  and  the  comparison  scene. 

The  ideas  demonstrated  by  Heflinger  et  al.  (1966)  have  advanced  the  field  of 
interferometry  for  use  in  aerodynamic  research  in  several  ways.  Since  reconstruction  of 
the  interference  pattern  required  a single  reference  beam  and  a single  common  path, 
precision  optical  components  and  accurate  mechanical  alignment  were  not  required.  The 
only  requirement  was  for  stability  between  the  two  exposures.  This  represented  a 
significant  step  forward  in  the  prospect  of  using  interferometric  measurements  in  the 
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harsh  conditions  associated  with  wind  tunnels  and  other  aerodynamic  measurement 
environments.  In  addition  to  these  benefits,  the  differential  interference  pattern  obtained 
with  this  method  meant  that  observations  could  be  made  through  windows  of  poor 


(a)  first  exposure  (b)  second  exposure  (c)  reconstruction  of  interferogram 


optical  quality.  The  work  published  in  this  paper  demonstrated  the  feasibility  of  the 
concept  by  making  interferograms  of  a mach  1.7  and  mach  2.5  projectile  as  well  as  of  the 
temperature  distribution  inside  a lighted  incandescent  lamp  bulb. 
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Holographic  interferometry  in  conjunction  with  computed  tomography 

In  1 966  Rowley  proposed  using  holographic  interferograms  along  with  the 
principles  of  reconstruction  from  projections  to  determine  the  three  dimensional  index  of 
refraction  of  a test  region.  As  discussed  above,  Matulka  and  Collins  (1971)  were  the  first 
to  demonstrate  the  use  of  holographic  interferometric  tomography  (HIT)  to  reconstruct 
the  density  field  around  a supersonic  free  air  jet.  Since  that  time  holographic 
interferometry  has  been  suecessfully  used  in  numerous  aerodynamic  measurements  (Zien 
et  al.  1974;  Snyder&Hesselink  1984;  Vukicevic  et  al.  1989;  Doerr  1990;  Mao  et  al.  1991; 
Timmerman  & Watt  1997). 

Phase  shifting/double  plate  method 

There  have  been  several  advances  in  holographic  interferometry  since  its  inception 
which  are  worth  mentioning  here.  Many  of  these  advances  are  concerned  with  the 
retrieval  of  phase  from  the  interference  fringe  pattern  and  apply  to  all  forms  of 
interferometry.  As  discussed  above  there  was  originally  no  attempt  to  control  the 
absolute  phase  of  the  reference  beam  between  exposures  or  at  reconstruction.  This  means 
that  the  interferogram  represents  a fixed  phase  difference  between  the  test  and 
comparison  beams,  which  cannot  be  changed.  In  1974  Abramson  proposed  a new 
method  called  Sandwich  Hologram  Interferometry  in  which  two  separate  holograms  were 
made  of  the  test  and  comparison  scenes  which  were  then  sandwiched  together  for 
reconstruction.  By  mechanically  altering  the  alignment  between  the  two  holograms 
during  reconstruction,  additional  reference  fringes  and  or  phase  shifts  can  be  introduced. 
In  a similar  technique  Dandliker  et  al.  (1976)  and  others  have  used  separate  reference 
beams  for  each  exposure  of  the  single  hologram  double  exposure  technique.  The  relative 
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phase  between  the  reference  beams  can  be  varied  during  reconstruction  to  achieve 
different  interferogram  phases.  The  relative  tilts  between  the  reference  beams  could  also 
be  varied  to  generate  reference  fringes.  In  both  of  these  techniques  the  desire  to  change 
phase  between  interferograms  stems  from  the  need  to  solve  for  the  phase  values 
associated  with  the  measured  intensity  values  of  the  interference  pattern.  Multiple 
interference  patterns  generated  with  different  relative  phases  will  result  in  as  many 
equations  relating  phase  to  the  interference  intensity  pattern.  By  measuring  3-4  such 
patterns,  all  ambiguities  can  be  eliminated  allowing  for  the  solution  of  the  actual  phase 
distribution. 

Fourier  transform  fringe  analysis 

Another  advance  came  about  in  1982  when  Takeda  et  al.  introduced  the  Fourier 
transform  method  of  phase  retrieval.  In  this  method,  carrier  fringes  are  added  to  the 
interference  pattern  and  the  principles  of  Fourier  analysis  are  used  to  separate  out  and 
solve  for  the  phase  distribution.  The  addition  of  carrier  fringes,  which  a simply  achieved 
by  tilting  of  the  reference  beam  between  exposures,  is  required  to  separate  the  modulation 
spectrum  from  the  DC  component.  The  method  is  straightforward  and  easily  applied 
using  fast  Fourier  transform  methods,  however  the  addition  of  the  carrier  fringes  limits 
the  resolution  of  the  reconstructed  phase. 

Limited  data  reconstruction 

In  many  measurement  applications  the  region  of  interest  contains  an  obstruction 
which  blocks  a portion  of  the  projection  measurement  and  in  some  cases  it  is  impossible 
to  obtain  projections  over  the  full  180  degrees  of  viewing  angle  required  for  ideal 
reconstruction  from  projections.  A number  of  authors  have  done  a significant  amount  of 


research  in  this  area  and  have  developed  methods  for  performing  reconstructions  with 
limited  data  (Sato  et  al.  1981;  Cha  & Sun  1989;  Cha  & Sun  1990). 

Weaknesses  of  HIT  techniques 
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Holographic  interferometry  has  proven  to  be  a valuable  tool  in  aerodynamics 
research  over  the  past  3 decades,  however  it  does  have  several  shortcomings  which  have 
prevented  its  wide  scale  use.  The  photographic  film  used  for  recording  the  holograms 
requires  processing  complexity  and  delay  and  is  one  of  the  main  vulnerabilities  of  this 
technique.  A great  deal  of  care  must  be  taken  in  exposing  the  film  in  order  to  avoid  non- 
linearities  and  in  maintaining  the  stability  of  the  holograms  over  time.  Although  it  made 
significant  strides  in  overcoming  the  alignment  and  precision  optics  requirements  of 
conventional  interferometry,  holographic  interferometry  setups  can  still  be  quite 
cumbersome  to  implement.  Typically  a minimum  of  9-10  different  projection  view 
angles  are  required  for  accurate  reconstruction.  This  represents  9 test  beams,  9 reference 
beams,  and  9 holographic  film  plates  routed  through  and  surrounding  a common  test 
region.  In  many  applications  this  kind  of  arrangement  becomes  prohibitively  complex. 
Because  of  these  weaknesses  there  is  an  interest  in  finding  other  methods  of  generating 
the  multidirectional  projection  data. 

Common  Path  Interferometry 

One  of  the  primary  drawbacks  of  both  the  direct  and  holographic  interferometry 
techniques  discussed  above  is  the  fact  that  a separate  reference  beam  is  required  in 
addition  to  the  object  beam  for  each  projection  angle.  In  many  test  environments  the 
difficulty  with  the  routing,  alignment,  and  vibration  control  for  many  different  beams  is 
prohibitive  and  these  techniques  become  impractical.  Phase  contrast  visualization  is  a 
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technique  originally  employed  in  the  viewing  of  phase  objects  in  the  field  of  microscopy 
(Zemike  1942;  Smartt  & Steel  1975).  It  has  since  been  extended  in  the  form  of  common 
path  interferometry  and  applied  to  flow  visualization  (Taylor  1980;  Anderson  & Taylor 
1982;  Paul  1982;  Bachalo  & Houser  1985;  Anderson  & Lewis  1985;  Loomis  et  al.  1991). 
There  are  many  variations  on  this  technique,  but  the  underlying  idea  is  to  eliminate  the 
need  for  a separate  reference  wave.  Common  path  interferometers  are  distinguishable  in 
that  they  use  a common  path  for  the  test  and  reference  beams  and  therefore  eliminate 
many  of  the  mechanical  alignment  and  vibration  problems  associated  with  dual  path 
interferometers  such  as  the  Mach-Zehnder.  By  using  a spatial  filter  in  the  Fourier  plane 
of  a simple  telescopic  optical  system,  which  ‘synthesizes’  a reference  beam  by  diffraction 
while  having  negligible  effect  on  the  object  beam,  an  interference  pattern  is  generated 
very  much  as  if  the  reference  beam  were  provided  separately.  This  method  works  well  for 
qualitative  visualization.  When  the  object  wave  has  low  spatial  frequency  content  the 
synthesized  reference  wave  can  be  assumed  to  be  planar  with  little  error.  However,  in  the 
completely  general  case  where  the  object  beam  can  have  high  spatial  frequency  content, 
the  spatial  filter  does  not  diffract  an  ideal  spherical  reference  wave  and  it  is  not  possible 
to  reconstruct  the  actual  phase  values  making  it  unsuitable  for  quantitative  analysis. 
Theory  of  oneration 

Point  diffraction  interferometers  operate  on  the  same  principles  derived  by  Zemike 
in  1942  for  use  in  phase  contrast  microscopy  and  are  similar  in  operation  to  Schlieren 
methods.  As  shown  in  Figure  2-9,  a spatial  filter  in  the  Fourier  plane  of  a telescopic 
optical  system  is  used  to  alter  the  phase  and/or  amplitude  of  a portion  of  the  spatial 
frequency  spectmm.  It  is  straightforward  to  show,  by  Fourier  analysis,  that  this  will 
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make  phase  objects  visible  in  the  output  image  plane.  Conceptually,  the  spatial  filter 
which  is  typically  the  size  of  the  system’s  diffraction  limited  spot  size,  diffracts  a 
spherical  wave  from  the  DC  focus  of  the  input  lens  while  passing  the  remainder  of  the 
spectrum  with  little  or  no  change.  The  diffracted  spherical  wave  can  be  thought  of  as  the 
reference  wave,  which  becomes  approximately  planar  following  the  second  lens  and 
interferes  with  the  unaltered  high  pass  portion  of  the  spectrum. 


One  common  PDI  technique  is  Dark  Central  Ground  which  employs  a spatial 
Fourier  filter  that  is  simply  a dark  spot  at  the  DC  or  central  location  with  a diameter  on 
the  order  of  the  diffraction  limited  spot  size.  There  are  numerous  other  techniques, 
which  differ  primarily  in  the  form  of  the  spatial  filter  used.  The  basic  Fourier  analysis  is 
as  follows  (Anderson  1995).  At  the  Fourier  plane  of  lens  one  the  amplitude  distribution 


can  be  written  as. 


27 


j'2-n. 


S(^,r)  = cjj s(x,y)e  ^ ^’'^dxdy  = 3 


circ 


+/ 

R 


[2.11] 


where  s(x,y)  is  the  signal  beam  after  passing  through  the  object 
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and  g(x,y)  represents  the  amplitude  modulation  of  the  beam  and  0(x,y)  represents  the 
phase  modulation  or  distortion  that  is  being  measured. 

Any  common  path  spatial  filter  function  can  be  written  in  the  following  form: 


H(^,r)  = a 


1 + be^*circ 


[2.13] 


where,  (^„,x<,)  is  the  center  of  the  spatial  filter  (i.e.,  the  location  of  the  diffracting  spot), 
and  Wf  is  its  radius.  The  parameters  a,  b,  and  <j)  determine  the  filter  type.  For  example, 
a=l,  b=l,  and  for  the  dark  central  ground  filter.  The  amplitude  and  phase 

distribution  at  the  filter  output  can  be  written  as 
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and  following  a Fourier  transform  by  the  second  lens  the  output  is 


o(x,y)  = a 


s{-x,-y)  + be^^3 


S{^,y)circ 


[2.15] 


The  first  term  in  equation  [2.15]  is  a reversed  version  of  the  original  object  wave 
while  the  second  term  represents  the  reference  wave.  The  superposition  of  the  two  terms 
in  equation  [2.15]  creates  an  interference  pattern  of  intensity  in  much  the  same  way  that  a 
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conventional  interferometer  does.  In  order  to  be  able  to  quantitatively  evaluate  the 
interference  pattern,  the  amplitude  and  phase  distribution  of  the  reference  wave  must  be 
known.  Recalling  that  a,  b,  and  ^ are  constants  set  by  the  form  of  the  spatial  filter,  the 
general  form  of  the  reference  wave  amplitude  and  phase  is  determined  by 


r{x,y)  = 3- 


S{^,r)circ 


or 


[2.16] 


r(x,y)  = 5(-.x,->’)0  3< 
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[2.17] 


Reference  wave  assumptions  and  dependence  on  object  field 

In  most  cases  Wf  is  very  small  as  compared  to  S(^,y)  so  the  second  term  in  the 
convolution  of  equation  [2.17]  is  very  flat  and  broad  which  has  a smoothing  effect  on 
5(-x,-y)  generating  a nearly  constant  amplitude  reference  wave.  If  the  spatial 
frequency  spectrum  is  smooth  so  as  to  be  nearly  constant  over  the  spatial  filter 

spot  located  at  then  the  phase  resulting  from  equation  [2.17]  will  also  be  nearly 

constant.  If  these  conditions  are  met  then  the  reference  wave  can  be  assumed  to  be  planar 
with  little  error.  The  reliance  on  these  conditions  is  the  primary  weakness  of  these 
methods  for  use  in  quantitative  analysis.  Since  the  object  field  s(x,y)  is  the  unknown 
subject  of  the  measurement,  the  assumptions  required  above  often  cannot  be  made.  In 
spite  of  these  shortcomings,  these  methods  have  been  used  often  for  flow  visualization 


measurements. 
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History  of  common  path  interferometry  in  fluid  flow  measurement 

One  application  in  which  the  assumptions  required  of  the  object  waye  often  hold  is 
in  the  measurement  of  aberrations  in  optical  components  and  systems.  The  common  path 
interferometer  is  especially  well  suited  for  these  measurements  since  it  is  readily  inserted 
in  the  optical  system  under  test.  Smartt  and  Steel  (1975)  and  Smartt  (1979)  reported  on 
the  use  of  point  diffraction  interferometry  in  measurements  on  astronomical  telescopes 
under  operating  conditions  with  a star  as  the  light  source  as  well  as  with  a laser  source. 
These  measurements  made  yisible  the  aberrations  and  system  misalignments  as  well  as 
perturbations  due  to  atmospheric  distortion.  Taylor  (1980)  reyiewed  the  principles  of 
operation  for  the  method  of  phase  contrast  as  applied  to  gas  flow  yisualization  and  made 
measurements  of  a gas  jet  flow.  He  recognized  the  shortcomings  of  the  required 
assumptions  and  pointed  out  that  the  method  is  more  suited  to  measurements  where  the 
density  yariations  are  small  resulting  in  sub-wayelength  optical  path  differences. 

Anderson  and  Taylor  (1982)  looked  at  methods  of  enhancing  the  fringe  pattern  accuracy 
by  tapering  the  spatial  filter  spot  at  the  edges  by  simulation  and  presented  some 
measurements  of  a gas  jet  flow.  Anderson  and  Milton  (1989)  reported  on  a large  aperture 
inexpensiye  interferometer  for  routine  flow  measurements  which  was  based  on  a dark 
central  ground  filter.  They  used  the  laser  source  focussed  onto  a photographic  plate  in 
the  Fourier  plane  to  make  their  dark  central  ground  filters  and  published  interferograms 
of  the  temperature  distribution  around  a candle  flame.  Loomis  et  al.  (1991)  reported  on 
the  application  of  dark  central  ground  method  as  published  by  Anderson  and  Milton  for 
use  in  supersonic  wind  tunnel  testing.  Their  paper  demonstrated  shock-shock  interactions 
as  well  as  general  flow  around  obstructions  in  a wind  tunnel.  Brock  et  al.  (1991) 
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examined  point  diffraction  interferometry  for  use  in  wind  tunnel  testing  and  presented 
some  measurements  taken  of  subsonic  flow  around  an  airfoil.  In  this  paper,  some 
advantages  of  PDI  over  DCGI  were  pointed  out  and  the  effects  of  pinhole  size  and  flow 
conditions  were  examined  qualitatively  in  terms  of  their  effect  on  fringe  visibility. 
Medecki  et  al.  (1996)  used  a point  diffraction  interferometer  in  a novel  arrangement  with 
a moveable  grating  to  produce  multiple  fringe  patterns  with  variable  phase  shifts.  This 
application  was  in  the  measurement  of  aberrations  of  optics  components  and  the  use  of 
phase  stepping  allowed  them  to  make  highly  accurate  measurements. 

Weaknesses  of  common  path  interferometry 

The  primary  limitation  of  the  common  path  methods  is  the  fact  that  the  intensity 
values  of  the  common  path  fringe  pattern  are  dependent  upon  the  object  of  the 
measurement  thus  requiring  benign  object  fields  in  order  to  fulfill  the  assumptions 
required  for  the  above  reference  wave  approximations.  This  limitation  makes  common 
path  interferometers  a poor  choice  where  accurate  quantitative  data  is  required  for  general 
aerodynamic  flow  fields,  as  in  the  case  of  tomographic  reconstruction. 

Wavefront  Sensors 

Wavefront  sensors  are  different  from  the  above  methods  in  that  they  effectively 
sample  the  gradient,  or  wavefront  slope,  of  a beam.  If  absolute  phase  is  desired,  then 
these  gradients  must  be  integrated  across  the  aperture.  This  type  of  sensor  has  found  a 
niche  in  adaptive  optical  systems  for  atmospheric  aberration  compensation  in 
astronomical  telescopes.  They  have  been  extensively  used  in  this  field  (Hudgin  1977; 
Furber  & Jordan  1997).  In  these  applications,  the  wavefront  sensor  measurements  are 
used  in  a closed  loop  control  system  along  with  a deformable  mirror,  which  is  adjusted  to 
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compensate  for  atmospherically  induced  fluctuations  in  the  wavefront  phase.  The  two 
primary  wavefront  sensors  which  have  been  used  are  the  lateral  shearing  interferometer 
and  the  Shack-Hartmann  sensor.  The  two  have  been  shown  to  be  equivalent  in  operation 
and  measurement  accuracy  by  Pennington  et  al.  (1994)  and  others.  Due  to  recent 
advances  in  the  ability  to  fabricate  Hartmann  lens  arrays,  the  Hartmann  sensor  is 
becoming  more  widely  used  (Shiono  et  al.  1987;  Borrelli  & Morse  1988;  Popovic  et  al. 
1988;  Jahns&Walker  1990;  Kwo  et  all 991;  Artzner  1992;  Brown  1994).  It  requires  less 
complexity  in  its  optical  setup  and  for  these  reasons,  will  be  considered  in  lieu  of  the 
lateral  shearing  interferometer  for  this  work. 

Theory  of  operation 

In  the  Hartmann  sensor,  the  object  wave  is  broken  into  a spatial  array  of  samples  or 
sub-apertures  by  an  array  of  lenslets  as  shown  in  Figure  2-10.  Each  sub-aperture  is 
brought  to  focus  on  a detector  array  positioned  in  the  back  focal  plane  of  the  lenslets. 
Since  the  lateral  position  of  the  focussed  spots  depends  upon  the  local  tilt  of  the  incoming 
wave,  the  array  of  spot  displacement  measurements  represents  the  wavefront  gradient 
sampled  at  the  locations  of  the  lens  array  sub-apertures. 

Applying  simple  geometric  optics  to  a single  lens  of  the  array  reveals  that  the  spot 
displacement,  , is  related  to  the  local  tilt  angle,  6»,. , and  the  focal  length  of  the  lenslets, 
/,  according  to  the  following  equation: 

- /tan(6»,)  . [2.18] 

Since  the  displacement  is  linearly  proportional  to  the  focal  length  of  the  individual 
lenses  making  up  the  array, /is  an  important  parameter  determining  the  sensitivity  of  the 
system.  Often  the  small  angle  approximation  is  used  and  the  equation  is  written  in  the 
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following  form: 


[2.19] 


As  seen  from  the  above  description,  the  resolution  of  the  Hartmann  sensor  is 
determined  by  the  sub-aperture  lenslet  spacing,  and  the  accuracy  is  determined  primarily 
by  the  ability  to  measure  the  diffraction  spot  centroid  positions  on  the  detector  array  in 
the  presence  of  noise.  The  task  of  designing  a Hartmarm  sensor  for  quantitative  flow 
measurements  can  be  broken  down  into  3 parts  consisting  of  lens  array  design,  spot 
centroid  measurement,  and  reconstruction  of  wavefront  phase  using  the  measured 
gradient  data.  There  has  been  a significant  amount  of  research  done  in  all  these  areas 
independently,  which  has  been  primarily  steered  toward  adaptive  optics  applications  for 
astronomical  telescope  compensation. 
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The  lens  array  samples  the  incoming  wave,  which  may  or  may  not  be  magnified 
from  its  size  in  the  flow  test  region.  Intermediate  magnification  or  demagnification  could 
be  used  to  match  resolution  requirements  to  available  lens  array  dimensions.  Expected 
wavefront  gradients  are  used  to  determine  the  required  focal  length  for  the  required 
sensitivity  and  accuracy.  Diffractive  as  well  as  refractive  lens  arrays  have  been  studied 
for  use  in  Hartmann  sensors  and  numerous  processes  from  photolithography  to  photo- 
thermal  techniques  have  been  used  to  produce  them.  The  general  requirement  is  to  be 
able  to  produce  uniform  arrays  with  small  sub-apertures,  high  fill-factor,  and  typically 
high  f-number.  Arrays  with  sub-aperture  diameters  in  the  20  micron  to  4 millimeter 
range,  fill  factors  approaching  100%,  and  f-numbers  from  near  unity  to  several  hundred 
are  commercially  available  from  a number  of  manufacturers. 

Spot  deflection  measurement 

The  design  of  the  lens  array  cannot  be  separated  totally  from  the  problem  of  spot 
centroid  measurement  on  the  detector  array.  The  lens  array  must  be  matched  to  the 
detector  array  geometry  for  the  desired  centroid  measurement  scheme.  Quadrant  centroid 
measurements,  which  require  four  detector  pixels  per  subaperture  have  been  used,  as  well 
as  multiple  pixel  schemes  using  many  detector  pixels  per  subaperture  (McMackin  et  al. 
1994;  Masson  et  al.  1995;  Neal  et  al.  1995;  McMackin  et  al.  1995;  Chen  et  al.  1996; 
Sapayo  & Truman  1997).  Although  the  research  was  not  done  with  a Hartmann  sensor  in 
mind,  there  is  a broad  base  of  information  available  concerning  the  accuracy  of  centroid 
measurement  using  a CCD  array.  Tyler  and  Fried  (1982)  developed  the  image  position 
error  associated  with  a quadrant  detector  as  a function  of  the  signal  to  noise  voltage  ratio. 
Winick  (1986)  determined  the  Cramer-Rao  lower  bounds  for  this  measurement  using  one 
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or  two  dimensional  detector  arrays.  Down  (1992)  has  analyzed  the  image  position  error 
for  focal  plane  arrays  by  assuming  poisson  statistics  for  object  and  background  shot  noise 
and  detector  dark  noise. 

Wavefront  phase  reconstruction 

The  problem  of  reconstructing  the  wavefront  phase  from  a set  of  wavefront  sensor 
measurements  has  been  studied  extensively  for  the  case  of  atmospheric  distortion  and 
optimum  solutions  based  on  various  techniques  have  been  derived.  Recall  that  the 
measured  spot  deviations  represent  the  localized  wavefront  gradient  averaged  across  the 
sub-aperture  associated  with  each  spot.  The  wavefront  reconstruction  becomes  one  of 
fitting  this  data  by  means  of  a least  squares  or  some  other  criteria  to  an  estimate  of 
wavefront  phase.  The  techniques  used  fall  into  one  of  two  classes.  The  first  class 
consists  of  zonal  methods  in  which  the  estimate  is  iteratively  matched  to  the  measured 
data  for  each  spatial  zone.  Zonal  techniques  were  used  by  Hudgin  (1977)  as  well  as  by 
Fried  (1977)  in  deriving  an  optimal  estimate  of  the  wavefront  phase  from  a set  of 
wavefront  sensor  measurements  defined  over  a square  array  and  using  a least  squares 
solution.  Both  authors  showed  that  the  mean  square  wavefront  error  due  to  measurement 
noise  depends  on  the  number  of  measurement  points,  N,  according  to  A + ln(N).  Noll 
(1978)  derived  the  same  result  analytically.  Hunt  (1979)  formulated  the  same  problem 
using  matrix  methods  and  developed  a more  efficient  reconstruction  scheme  in  terms  of 
its  rate  of  convergence.  Cubalchini  (1979)  examined  the  reconstruction  problem  from  a 
modal  point  of  view,  which  is  the  second  class  of  techniques  used.  In  this  method  the 
estimate  is  based  on  a set  of  orthogonal  basis  functions,  which  can  be  represented  by  a set 
of  coefficients.  As  with  most  of  the  research  conducted,  Zemike  polynomials  were  used 
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as  the  set  of  basis  functions  because  they  are  well  matched  to  atmospheric  aberrations  and 
astronomical  telescope  applications.  Southwell  (1980)  and  Herrmarm  (1980)  examined 
the  reconstruction  problem  from  both  points  of  view  and  made  comparisons  between 
them.  Southwell  (1980)  showed  that  zonal  methods  are  superior  in  terms  of  noise 
propagation  as  well  as  being  more  computationally  efficient  and  that  the  superiority  is 
enhanced  with  an  increasing  number  of  measurements. 

History  of  wavefront  sensor  use  in  flow  measurement 

As  stated  above,  most  of  the  work  done  with  Hartmann  sensors  has  been  in  terms  of 
adaptive  optics  systems  for  atmospheric  compensation.  In  this  field  the  characteristics  of 
the  wavefront  aberrations  are  well  known  and  are  somewhat  benign  as  compared  to  many 
aerodynamic  flow  measurement  problems.  However,  there  has  been  some  work  done 
related  to  general  aerodynamic  flow  visualization.  Johnson  (1993)  performed 
simulations  of  the  tomographic  reconstruction  process  using  simulated  Hartmann  sensor 
wavefront  data  for  a gaussian  density  profile.  The  simulations  were  run  while  examining 
the  effects  of  projection  sampling  and  number  of  projections  used.  For  the  gaussian 
profile  used,  sub-aperture  spacing  of  less  than  two  sub-apertures  per  2cr  width  provided 
insignificant  error  reduction.  The  effects  of  limited  projection  data  were  also  examined 
and  the  simulation  showed  that  more  than  18  projection  angles  spaced  over  180  degrees 
were  ineffective  at  decreasing  reconstruction  errors.  Pederson  (1993)  performed  an 
experimental  verification  of  Johnson’s  work  by  using  a lateral  shearing  interferometer  to 
execute  tomographic  reconstruction  of  turbulence  around  a jet  nozzle.  The  work  of 
Johnson  and  Pederson  was  later  summarized  by  Roggeman  et  al.  (1995).  McMackin  et  al. 
(1994),  and  Masson  et  al.  (1995)  described  a system  developed  by  Phillips  Labs  using  a 
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one-dimensional  Hartmann  sensor  and  tomographic  analysis  for  measurements  of  the 
flow  field  of  a heated  air  jet.  This  system  used  a one  inch  long,  40-element  cylindrical 
lens  array  and  a 2048  pixel  linear  CCD  array  to  take  one  dimensional  ‘slices’  of  the  flow 
field  at  a several  kilohertz  frame  rate.  Since  the  flow  velocity  was  slow  (5  m/s)  these 
slices  could  be  stacked  to  generate  a two  dimensional  picture  of  the  flow  passing  the 
sensor  location.  Neal  et  al.  (1995)  and  McMackin  et  al.  (1995)  reported  on  an 
enhancement  to  the  Phillips  Lab  system  whereby  binary  optics  lens  arrays  were  used 
having  40,  75,  and  150  lenses  over  the  one  inch  length.  The  same  2048  pixel  liner  CCD 
array  was  used.  Chen  et  al.  (1996)  described  further  enhancements  to  the  Phillips  lab 
system  (POTTOS,  Fast  Optical  Tomography  of  Turbulent  Organized  Structures)  in  which 
eight  source-sensor  pairs  were  used  to  generate  the  multiple  view  angle  data  required  for 
tomographic  reconstruction.  Here,  64  element  lens  arrays  and  2048  pixel  detector  arrays 
were  used  in  linear  configurations,  and  a frame  rate  of  5 kHz  provided  time  resolved  data 
of  the  flow  of  a heated  air  jet.  Tomographic  reconstruction  was  performed  using  the 
algebraic  reconstruction  technique  (ART),  and  results  are  presented.  Sapayo  (1997) 
reported  on  measurements  concerning  the  development  of  axisymmetric  and  helical 
modes  in  heated  air  jets  using  the  Phillips  Lab  POTTOS  system. 

Weaknesses  of  wavefront  sensor  methods 

The  wavefront  sensor  methods  add  an  additional  level  of  complexity  to  the 
reconstruction  process  since  the  projection  data  (wavefront  phase)  must  be  reconstructed 
from  the  spot  deflection  measurements.  Care  must  therefore  be  taken  to  ensure  that  noise 
in  the  spot  deflection  measurements  does  not  significantly  degrade  the  reconstruction. 
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Design  Example 

Applying  the  Gladstone-Dale  relationship,  the  supersonic  conical  flow  model,  and 
the  Hartmann  sensor  fundamentals  to  the  design  example  should  provide  some  additional 
insight  into  the  principles  covered  in  this  chapter.  Specifically,  it  should  be  possible  at 
this  point  to  establish  some  guidelines  for  the  anticipated  dynamic  range  of  the 
measurement. 

Maximum  Wavefront  Slope 

Note  that  we  have  developed  a solution  for  the  integrated  phase  data  as  a function  of 
the  field  point  angle,  0.  We  have  used  numerical  integration  to  solve  for  the  absolute 
phase  shift  at  a normalized  value  of  Zo  such  that  rs=l  cm,  as  0 is  varied  uniformly  from  0c 
to  0s.  For  an  arbitrary  slice  location  Zq,  the  phase  is  a simple  linear  scaling  of  the 
normalized  value  because  the  index  is  dependent  only  on  0.  As  seen  in  Figure  2-4,  the 
wavefront  slope  increases  to  infinity  at  the  shock  boundary.  At  this  boundary  the  input 
beam  will  be  strongly  refracted  such  that  the  rays  passing  through  this  immediate  vicinity 
will  fall  outside  of  the  field  of  view  of  the  measurement.  As  with  all  other  flow 
visualization  methods,  the  region  adjacent  to  the  shock  boundary  cannot  be  directly 
measured.  Therefore  it  is  logical  to  step  away  from  this  boundary  when  calculating  the 
maximum  anticipated  wavefront  slope  to  be  used  as  a design  specification.  The  step  size 
used  for  0 in  the  above  calculations  was  0.1°.  For  the  50mm  square  test  region,  this 
angular  step  size  corresponds  to  a linear  step  size  of  50-90  micrometers,  which  is  on  the 
same  order  as  a single  lenslet  of  a Hartmann  lens  array.  Therefore,  we  will  use  the  phase 
points  at  0 = 0s-O.l°  and  0 = 0s-O.2°  for  the  calculation  of  the  maximum  anticipated  phase 
slope.  Recall  that  we  must  scale  the  normalized  phase  shift  values  by  the  actual  shock 
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radius  at  the  slice  plane  defined  by  z=Zq.  From  Figure  2-1,  it  is  clear  that  the  shock  radius 
for  an  arbitrary  Zq  can  be  written  as  follows: 

. [2.20] 

This  can  be  expressed  in  terms  of  the  radial  distance  from  the  vertex  of  the  cone  to  the 
field  point  and  the  field  point  angle  0 as  follows: 

r^=r^ian6^cos9  . [2.21] 

The  absolute  phase  shift  (j){r,6)  of  a ray  passing  through  a field  point  at  {r^,  0 ) will  be 

the  normalized  phase  shift  ^[9]  scaled  by  this  value  of 

<!>{r„,9)  = r^ian9^cos9^^{9)  . [2.22] 

The  linear  distance  in  the  0 direction  between  the  last  two  points  prior  to  the  shock, 
which  are  located  at  (ro,  0 ) and  (ro,  0 ) can  be  closely  approximated  for  small  angles  as 
follows: 


M = rSe^-9,)  = r„^9  . [2.23] 

Now  the  phase  gradient  between  these  two  boundary  point  can  be  found  as  follows: 

1^(1,  tan  9^  [^(^2 ) cos(6»2  ) - ^(6>, ) cos(6>, )] 


r = 


A/ 


^9 


[2.24] 


Note  that  the  gradient  is  independent  of  Vo,  In  order  to  relate  the  phase  gradient  to  a 
wavefront  tilt,  we  must  simply  transform  the  phase  shift  to  its  corresponding  distance  in 
the  undisturbed  medium,  which  is  accomplished  by  a simple  division  by  the  wave 
number,  ko.  Therefore  the  wavefront  tilt  angle  can  be  written  as  follows: 


© = tan' 


J 


- tan' 


tan  9^  \^{9^ ) cos{9^ ) - ^(6>, ) cos(6>, ) j 


k^9 


[2.25] 
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Given  the  conditions  of  the  design  example  (Mach  1 .5,  0^  -60°,  A = 532nm) , the 
numerical  solution  results  in  normalized  phase  values  for  the  last  three  points  prior  to  the 
shock  boundary  of -131.1°,  -185.7°,  and  -226.6°.  These  points  are  separated  in  angle  by 
A0  =0.1°.  The  resulting  maximum  wavefront  tilt  is  therefore 


If  we  were  to  step  back  another  0.1°  from  the  shock,  using  the  field  points  at  6 =59.8°  and 
59.7°,  the  resulting  maximum  wavefront  tilt  would  be  0.18°. 

Verification  of  Small  Refraction  Assumption 

Given  this  initial  estimate  of  the  maximum  anticipated  wavefront  tilt,  it  is  imperative  that 
we  examine  the  fundamental  'small  refraction'  assumption.  This  assumption  simply 
requires  that  the  rays  passing  through  this  region  of  maximum  index  gradient  will  not  be 
refracted  beyond  the  range  of  the  measurement  capabilities.  This  requires  that  the 
maximum  deviation  of  these  rays  , Ax,  be  smaller  than  a lenslet  aperture  when  imaged 
onto  the  lens  array  as  shown  in  Figure  2-11.  Considering  the  test  region  size  of  50mm 
and  a typical  off-the-shelf  lens  array  size  of  25mm,  it  is  anticipated  that  a magnification 
stage,  M,  of  approximately  0.5  will  be  required.  Since  the  maximum  index  gradient, 
which  causes  the  maximum  wavefront  tilt,  occurs  at  the  shock  boundary,  the  rays  passing 
through  this  region  will  experience  all  of  their  refraction  at  this  point,  which  lies  near  the 
center  of  the  test  region.  These  rays  will  then  travel  a distance  of  //2=25mm,  at  the 
maximum  wavefront  tilt  angle,  to  the  output  of  the  test  region.  The  deviation  at  the 
image  plane  can  then  be  easily  approximated  as  follows: 


0 


max 


= 0.23°  . [2.26] 
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AXf  = MAx  = tan  = 0.5  x 25mm  x tan(0.23°)  - - SOjum  . [2.27] 


Figure  2-11:  Geometry  for  small  refraction  assumption 


Under  the  above  conditions,  the  small  refraction  condition  requires  that  the  lenslet 
aperture  size  be  larger  than  50pm. 

Summary 

In  this  chapter  we  have  examined  a broad  range  of  flow  visualization  and 
quantitative  flow  measurement  techniques,  and  have  provided  some  insight  into  the 
motivation  for  the  use  of  the  wavefront  sensor  for  the  current  work.  In  summary,  the 
wavefront  sensor  provides  a means  for  quantitatively  measuring  the  projection  data 
required  for  a tomographic  reconstruction  of  the  three-dimensional  fluid  density  profile. 
Because  the  wavefront  sensor  does  not  require  a separate  reference  beam  as  the 
interferometric  methods  do,  its  use  simplifies  the  optical  complexity  and  alignment.  This 
is  especially  important  in  the  present  application  since  the  tomographic  reconstruction  of 
a high-speed  aerodynamic  event  will  require  multiple  simultaneous  projection 


41 


measurements.  Since  the  wavefront  sensor  measures  localized  wavefront  slope  as 
opposed  to  a direct  phase  measurement,  its  use  requires  an  additional  step  of 
reconstruction  to  obtain  the  phase  measurement.  A careful  analysis  of  the  errors 
introduced  in  this  process  must  be  performed.  Figure  2-12  is  a generalized  block  diagram 
showing  the  steps  required  for  the  use  of  the  wavefront  sensor  as  the  projection  data 


Figure  2-12:  Flow  diagram  of  measurement  process 


measurement  system  of  a tomographic  reconstruction.  Recall  that  the  final  steps 
involving  the  tomographic  reconstruction  will  not  be  performed  in  the  measurements 
since  only  a single  projection  measurement  system  is  being  built.  However,  the  system  is 


being  designed  and  built  under  the  assumption  that  a tomographic  reconstruction  would 
be  the  final  product.  Therefore  the  requirements  of  such  a system  will  be  considered 
carefully  in  the  design  of  this  projected  phase  measurement  system. 

The  supersonic  conical  flow  model,  which  is  being  used  as  the  design  example 
throughout  this  work  indicates  that  the  maximum  anticipated  wavefront  slope  will  be  on 
the  order  of  0.2°.  Due  to  the  abrupt  discontinuity  in  fluid  density  at  the  shock  boundary, 
the  region  adjacent  to  the  shock  will  not  produce  valid  measurements.  A review  of  the 
small  refraction  conditions  which  must  be  assumed  for  the  Shack-Hartmann  wavefront 
sensor  measurement  indicates  that  the  condition  is  met  for  wavefront  tilt  angles  which  are 
less  than  or  equal  to  0.2°  as  long  as  the  lens  array  sub-aperture  size  is  greater  than  about 
50pm. 


CHAPTER  3 

TOMOGRAPHIC  RECONSTRUCTION 
OF  THREE-DIMENSIONAL  DENSITY 

Problem  Statement 

In  this  chapter  we  will  examine  the  tomographic  reconstruction  process  which 
would  be  required  to  perform  the  three-dimensional  density  reconstruction.  Since  we  are 
collecting  only  a single  projection  angle  of  measurement  data,  the  tomographic 
reconstruction  will  not  be  performed  in  this  work.  However,  it  is  important  to  understand 
the  process  in  terms  of  reconstruction  accuracy  and  noise  performance  to  establish  a set 
of  requirements  for  the  projection  data.  With  this  in  mind,  we  will  examine  the  various 
tomographic  reconstruction  algorithms,  which  are  suited  to  the  type  of  measurements  we 
will  be  performing.  We  will  then  look  at  the  effects  of  discrete  sampling,  such  as  by  a 
Shack-Hartmann  sensor,  and  then  perform  a noise  analysis  in  order  to  understand  how 
projection  data  noise  affects  the  tomographic  reconstruction. 

Much  of  the  following  analysis  is  based  upon  the  convolutional  back-projection 
algorithm,  which  is  but  one  of  a number  of  tomographic  reconstruction  algorithms. 

Many  of  the  algorithms  used  in  this  type  of  work  are  iterative  and  as  a result,  a closed 
form  analysis  of  their  performance  is  not  possible.  Since  the  convolutional  back- 
projection  algorithm  is  one  of  the  fundamental  techniques,  the  results  obtained  from  its 
analysis  will  provide  some  insight  into  the  properties  of  the  other  methods.  Apart  from 
this  type  of  analysis,  there  is  empirical  data  available,  which  establishes  guidelines  for  the 
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requirements  of  the  projection  data,  and  these  results  will  be  presented  as  well.  In 
addition,  a representative  iterative  reconstruction  technique  known  as  algebraic 
reconstruction,  ART,  will  be  developed  in  Matlab.  It  can  then  be  used,  along  with  the 
design  example,  to  examine  the  expected  accuracy  and  noise  performance. 

Background 

Tomographic  Reconstruction  of  Three  Dimensional  Refractive  Index 

The  basic  problem  of  reconstructing  a three-dimensional  density  field  from  its 
projections  lies  in  the  inversion  of  the  line  integrated  optical  path  lengths.  For  a given 
viewing  angle, y , and  a given  cross  sectional  plane,  , the  measured  projection 

data,  («,'),  represents  the  integrated  optical  path  lengths  incurred  by  an  optical  beam 

traversing  the  medium  with  index  of  refraction  «(m,  , ) . Refer  to  Figure  3- 1 . 

00 

Py{u^')=  J«(M,'cosy -M2'siny,M,'siny + M2i'cosy)c/M2'  [3.1] 


Figure  3-1:  Projection  coordinate  system 
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The  projection  data  function,  p^{u^') , represents  the  integrated  optical  path  length  due 
to  a cross  sectional  slice  of  a plane  wave  traversing  a medium  with  unknown  index  of 
refraction.  To  reconstruct  the  two-dimensional  function,  «(m,  ,«2) , for  the  general  case 
where  no  symmetry  is  assumed,  a full  set  of  projection  data  is  required  by  varying  y 
over  180  degrees.  The  problem  of  inverting  equation  [3.1],  to  solve  for  n{u^,U2)  from  its 
projections  , Py{u^'),  is  known  as  computed  tomography  (CT)  or  computer  assisted 
tomography  (CAT).  It  is  most  well  known  for  its  application  to  the  field  of  medical 
radiology  for  which  the  Nobel  Prize  was  awarded  in  1979.  (Kak  & Slaney  1988) 

In  our  unique  application,  the  projection  data  represents  the  integrated  optical  path 
length  due  to  an  optical  ‘probing  beam’  as  opposed  to  representing  the  integrated 
absorption  due  to  x-ray  probing  as  used  in  the  medical  application.  In  theory,  the 
problem  of  inverting  equation  [3.1]  can  be  solved  exactly  using  the  Radon  transform. 
However  this  requires  that  the  projection  function,  p^  («,'),  be  defined  in  functional  form 

for  all  projection  angles.  Since Py{u^')  is  typically  sampled  discretely  in  the  m/ 
dimension  and  is  available  for  a finite  number  of  projection  angles,  the  analytical  solution 
cannot  be  used  and  various  numerical  techniques  have  been  developed.  One  exception  to 
this  limitation  is  in  the  case  of  axisymmetric  distributions  where  a single  projection 
represents  all  other  projection  angles  and  the  Radon  Transform  reduces  to  the  Abel 
Transform,  which  can  be  applied  to  yield  a direct  solution.  (Kak  & Slaney  1988) 

In  the  first  published  demonstration  of  three-dimensional  density  field 
reconstruction,  Matulka  and  Collins  ,(1971)  used  holographic  interferometry  and  series 
expansion  methods  to  reconstruct  the  density  field  around  a supersonic  free  air  jet.  The 
tests  demonstrated  the  ability  to  perform  this  measurement  with  an  axially  symmetric  as 
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well  as  a tilted  jet  asymmetric  case.  Their  use  of  holographic  interferometry  built  upon 
the  work  of  Heflinger  et  al.,(1966).  The  three  dimensional  reconstruction  methods  were 
based  on  an  extensive  foundation  of  work  done  in  the  field  of  plasma  physics  (Bockasten, 


Maldanado  et  ah,  (1966)  for  measuring  plasma  emissivity. 

The  central  slice  theorem 

Many  of  the  techniques  used  to  perform  the  reconstruction  from  projections  are 
based  on  the  ‘central  slice  theorem’.  This  theorem  states  that  the  one-dimensional 
Fourier  transform  of  the  projection  function,  Py{u^') , is  the  two-dimensional  Fourier 
transform  of  the  original  unknown  function,  n{u^,U2 ) A^(Q,  ^ evaluated  along  a 

line  which  passes  through  the  origin  at  an  angle  y from  the  transformed  coordinate 
Q,  (Kak&Slaney,  1988).  Using  the  coordinate  transform. 


1961  Barr  1962)  and  the  specific  numerical  method  was  taken  from  the  work  done  by 


M,  = M|'cosy  -Mj'siny 
= M|'siny  -i-Mj'cosy 


[3.2] 


Py{u^')  can  be  written  as 


—CO 


Taking  the  Fourier  transform  of  p^  (m,  ')  yields 


—CO 
00  CO 


[3.4] 


In  terms  of  the  un-rotated  coordinates. 
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00  00 


or 


[3.5] 


S^{o))  = N{coQ,osy,Q)s\nY)  . 


Equation  [3.5]  is  a mathematical  statement  of  the  central  slice  theorem. 

Using  this  theorem,  the  Radon  inversion  formula  can  be  derived.  The  unknown 
function  «(m,,M2)  can  be  found  by  inverse  Fourier  transforming  A^(Q,,Q2). 


Nvquist  sampling  requirements 

In  practical  measurement  applications,  the  projection  data  is  discrete  in  both  the 
spatial  and  angular  sense.  Since  the  projection  samples  are  linearly  spaced  by  an  equal 
increment  of  a , then  to  satisfy  Nyquist  sampling  requirements,  the  projection  data  should 
be  bandlimited  to  < ^ . Assuming  a total  spatial  extent  for  the  projection  data  of  d , the 
radial  Fourier  transform  samples  will  be  separated  by  . Assuming  that  M evenly 
spaced  projection  angles  over  180  degrees  are  used,  the  radial  spokes  in  the  Fourier  plane 
will  be  separated  by  radians,  as  shown  in  Figure  3-2.  In  building  up  the  two- 
dimensional  Fourier  space  using  M evenly  spaced  projection  angles,  the  outermost  radial 
Fourier  domain  samples  will  be  spaced  by. 


[3.6] 


In  terms  of  polar  coordinates 


[3.7] 


[3.8] 


By  equating  this  to  the  radial  spacing,  we  can  generate  a requirement  for  the  number  of 
projection  angles. 
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In  ^ 

= — M> 

Mad  2a 


[3.9] 


Assuming  that  we  have  2N+1  samples  evenly  spaced  over  the  full  extent  d , then  a can 
be  written  in  terms  of  d and  substituted  into  equation  [3.9]  to  obtain, 

M>Nn  [3.10] 


The  central  slice  theorem  and  Fourier  transform  domain  method  discussed  so  far 
provide  the  background  principles  upon  which  the  tomographic  reconstruction  process  is 
based  and  provide  an  intuitive  view  from  which  the  valuable  rules  of  thumb  derived 
above  can  be  obtained.  However  the  variable  frequency  resolution  obtained  from  the 
polar  format  of  the  projection  data,  and  the  interpolation  required  to  fit  this  data  to 
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cartesian  coordinates  prior  to  the  final  inverse  Fourier  transform  is  a fundamental 
problem  of  this  method  which  has  prevented  its  widespread  use.  There  are  numerous 
other  methods  which  are  more  often  used  in  practice,  and  which  can  often  provide  better 
results  in  terms  of  reconstruction  accuracy  and  computational  efficiency. 

Convolutional  back-projection 

One  of  the  most  commonly  used  CT  algorithms  is  the  convolutional 
backprojection  algorithm  also  known  as  filtered  backprojection,  (Kak&Slaney,  1988).  It 
is  used  extensively  in  medical  CT  applications  and  any  time  that  full  sets  of  projection 
data  are  available.  Rewriting  equation  [3.7]  in  a different  form, 

1 f 

= — g (m,  cosy  + «2  s\ny)dy  [3.11] 

2k  { 

where,  gy{t)  = Py{t)®  k{t)  [3.12] 

and  k(t)  is  the  inverse  Fourier  transform  of  |u;|  yields  a new  interpretation.  Here,  g^  (t) 
is  the  result  of  filtering  the  1-D  projection  data,  (t)  , with  a filter  whose  frequency 
response  is  \co\ . In  discrete  form,  equation  [3.11]  can  be  approximated  as 

M-\ 

sM^cosy + u^smy)  . [3.13] 

i=0 

Effectively,  the  filtered  projections  are  smeared  back  across  the  2-D  space  along 
the  direction  from  which  they  were  projected  and  averaged  with  the  filtered  projections 
taken  from  all  other  angles.  The  term  back-projection  comes  from  the  fact  that  the 
function  g,(Mj  cosy  + U2  siny)  = g,(M,') , which  is  uniform  in  , is  one-dimensional  and 
is  used  to  obtain  the  two  dimensional  function.  Although  the  ideal  frequency  response  of 
the  filter  k(t)  is  \o)\ , the  projections  are  assumed  to  be  band-limited  so  that  the  high 
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frequency  behavior,  which  would  amplify  high  frequency  noise,  is  neither  necessary  nor 
beneficial.  For  this  reason  the  filter  is  normally  chosen  to  be  approximately  linear  out  to 
some  cut-off  frequency  after  which  it  falls  to  zero.  The  cut-off  frequency  should  be  set  at 
no  more  than  one-half  of  the  sample  frequency  to  avoid  aliasing. 

Projection  data 

In  the  intended  application,  the  projection  data  will  be  measured  optically,  by 
passing  a coherent  plane  optical  beam  through  the  test  region  and  measuring  the 
accumulated  phase  shift.  The  phase  shift  will  be  sampled  over  a two-dimensional  grid, 
which  is  oriented  orthogonally  to  the  beam  direction,  with  a sample  spacing  of  a and 
having  an  odd  number  (2N+1)  of  points  in  both  dimensions.  The  data  will  be  separated 
into  one  dimensional  slices  of  length  2N+1,  each  of  which  make  up  a separate  set  of 
projection  data  to  be  used  in  the  tomographic  reconstruction.  The  reconstructed  slices 
will  then  be  combined  to  reconstruct  the  two-dimensional  reconstruction. 

The  phase  shift  experienced  by  a wave  passing  through  a media  of  index  n over  a 
length  of  travel  d,  can  be  expressed  in  terms  of  the  wave  number  ko=27i/A,  as  follows. 

A^  = k^nd  [3.14] 

If  the  index  is  a function  of  position  along  the  path  of  travel  /,  then  the  phase  shift  is  a 
function  of  the  integration  over  / of  the  index,  as  follows. 

Acf)  = k^^n{x)dx  [3-15] 

I 

Convolution  of  projection  data 

As  stated  in  equation  [3.13],  the  first  step  is  to  filter  the  projection  data  with  a 
filter  whose  frequency  response  is  \co\ . This  filter  is  normally  band-limited  as  discussed 
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in  section  2 to  one-half  the  sampling  frequency  l/2a.  A plot  of  the  frequency  response  is 
shown  in  Figure  3-3. 


- l/2a 

C 

l/2a 

Cm=l/2a 

Figure  3-3:  Backprojection  filter  frequency  response 


In  order  to  implement  the  convolution  operation,  we  need  to  analytically  obtain  the 
impulse  response  by  inverse  transforming  the  above  frequency  response.  This  is 
straightforward  and  results  in  the  following  continuous  impulse  response. 


h{x)  = 


smc 


smc 


2a 


[3.16] 


The  projection  data  is  discretely  sampled  at  x = na  where  n is  an  integer  ranging 
from  -N  to  N.  Substituting  for  x yields  the  discrete  impulse  response,  which  is  needed 
for  the  discrete  convolution  operation.  The  resulting  impulse  response  was  first 
described  by  Ramachandran  and  Lakshminaryanan  ,(1971),  and  will  thus  be  referred  to  as 
hRL(n). 


V(«) 


^ 1 ^ 


sinc(n)  - 


smc 


[3.17] 


Another  popular  version  of  the  filter  uses  a sinc(^2  C,  m)  taper  function  and  was 
developed  by  Shepp  and  Logan,  (1974).  Figure  3-4  shows  its  frequency  response. 
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_ l/7ta 
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^m=l/2a 

Figure  3-4  : Shepp-Logan  backprojection  filter  frequency  response 


The  impulse  response  for  the  Shepp-Logan  filter  is  as  follows. 


r 2 ^ 

1 

\7t  a J 

_4«'  -1_ 

[3.18] 


Backprojection 

Now  that  we  have  filtered  projection  data  for  M different  projection  angles,  the 
backprojection  function  as  described  by  equation  [3.13]  must  be  performed.  In  this 
process  we  must  sum  up  the  contributions  from  each  projection  angle  over  a fixed  grid  of 
reconstruction  points.  The  rotated  projection  data  will  not  line  up  with  the  un-rotated 
reconstruction  grid  so  a two  dimensional  interpolation  is  required.  The  rotated 
coordinates  of  the  projection  data  must  be  transformed  to  the  reconstruction  grid 
coordinates  using  the  simple  coordinate  rotation  transform  as  follows: 


xr  = xp  cosy  - yp  sin  y 
yr  = xpsmy  + ypcosy. 


[3.19] 
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The  one  dimensional  projection  data,  for  each  projection  angle,  is  first  back- 
projected  to  two  dimensions,  which  simply  involves  duplicating  the  projection  data 
vector  for  all  columns  of  a square  matrix.  This  data  is  then  rotated  according  to  equation 
[3.19]  and  the  points  are  interpolated  onto  the  un-rotated  reconstruction  grid  to  generate  a 

set  of  rotated  back-projections,  bakpro  ^ . These  steps  are  repeated  for  each  projection 
angle  and  accumulated  to  reconstruct  the  data.  The  reconstruction  will  be  of  the  following 
form. 


The  projection  data  as  given  in  equation  [3.12]  must  be  discretized  in  order  to 
carry  out  the  reconstruction  shown  by  equation  [3.13].  In  our  application,  this  sampling 
is  accomplished  by  the  Shack-Hartmann  sensor,  which  performs  an  average  wavefront  tilt 
measurement  over  each  lenslet  aperture.  It  is  important  to  consider  the  low  pass  filtering 
effects  of  this  sample  averaging  operation,  which  are  inherent  in  the  Shack-Hartmann 
sensor  measurement. 

Finite  Detector  Width  Effects 

Since  the  projection  data  is  not  ideally  sampled  by  delta  functions,  we  must 
examine  the  low  pass  filtering  effects  caused  by  the  finite  detector  width.  In  our 
application,  the  detector  width  is  set  by  the  size  of  the  individual  hartmaim  sensor 
lenslets.  Since  the  lenslet  effectively  averages  the  wavefront  tilt  across  its  aperture,  this 
averaging  operation  can  be  represented  by  a convolution  with  a uniform  detector  aperture 


[3.20] 


Sampling  Effects 


function  and  a division  by  the  width  of  the  aperture  function.  If  p(x)  is  the  continuous 
projection  data  and  a(x)  is  the  sampling  aperture  function,  where 
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a{x)  = 


1 I I ^ 

1 , X < — 

' ' 2 

0 , elsewhere 


[3.21] 


then  the  discrete  projection  data  can  be  written  as  follows: 


= 5{x  -na) 


W 


[3.22] 


Now  since, 


sm 


a{x)  « A{co)  = W- 


(OWy 


COWy 


[3.23] 


and. 


d{x)  = ^S{x  - na)  D{co)  = y'  s(o) 


Inn, 


[3.24] 


we  can  write. 


p„=3-’ 


P{(o)  A(ca) 
fV 


[3.25] 


As  seen  from  equation  [3.25],  the  projection  data  is  effectively  filtered  by  the  low  pass 
function. 


IF 


[3.26] 


Now  if  P(co)  is  bandlimited  to  ^ as  assumed,  and  the  sampling  window  has  a 100% 
duty  factor,  such  that  IV=a,  then  at  the  band  edges,  co  = ±y„ , the  low  pass  filter  of 
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equation  [3.26]  has  a value  of  y„ . This  is  the  form  used  in  the  Shepp-Logan  filter  and 
shown  in  Figure  3-4. 

The  effect  of  finite  sample  size  is  to  low  pass  filter  the  projection  data  as  described 
above.  The  degree  of  filtering  depends  upon  the  sampler  width.  At  a maximum  sampler 
duty  factor  of  1 00%,  the  sampling  implements  the  Shepp-Logan  filter  described  by 
equation  [3.18].  As  a result,  for  high  duty  factor  sampling,  the  convolution  step  of  the 
convolutional  back-projection  algorithm  need  only  implement  the  ideal  linear  frequency 
response  referred  to  above  as  the  Ramachandran-Lakshminarayanan  filter  to  obtain  the 
Shepp-Logan  reconstruction. 

Iterative  Reconstruction  Algorithms 

In  many  applications,  the  projection  data  is  limited  by  obscurations  within  the  test 
area  or  by  restrictions  on  the  available  projection  angles.  In  these  cases,  the 
reconstruction  can  be  performed  iteratively,  and  a number  of  algorithms  for  doing  so 
have  been  developed.  Most  iterative  reconstruction  algorithms  are  based  on  the 
principles  introduced  by  Gordon,  et  al.  (1970)  in  the  form  of  algebraic  reconstruction 
techniques  (ART).  The  ART  algorithms  are  more  intuitive  than  the  analytical  forms 
discussed  previously,  leading  to  a straightforward  although  computationally  intensive 
implementation.  The  basic  ART  algorithm  can  be  explained  as  follows. 

ART 

Let  the  reconstruction  space  be  divided  into  N non-overlapping  elements  and 
the  projection  data  will  also  be  divided  into  M non-overlapping  elements.  Let 
represent  the  projection  data  elements.  Let  represent  the  region  of  ^for  which  S\  is 
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the  projection.  Refer  to  Figure  3-5.  Now  if  we  let  r represent  a general  point  in  ^and 
rt(r ) is  the  unknown  index  function,  then  we  can  write  the  projection  as 

P.^k^\n{r)dr.  [3.27] 


Figure  3-5:  ART  geometry  and  space  representation 


The  intersection  of  with  the  i‘^  element  of  ^ can  be  defined  as  dP.  n which  allows 
us  to  write  the  projection  data  as 


\n{r)d7.  [3.28] 

'■=1 

In  performing  the  discrete  reconstruction,  we  will  be  assigning  a single  value  «,  to  the 
unknown  index  function  within  each  element  and  so  the  integral  in  [3.28]  may  be 
estimated  by  the  geometric  fraction  of  the  intersected  area  multiplied  by  the  discrete 


index  value  there,  such  that 
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I n{r)  dr  « n.  w,,  . 

dr 


[3.29] 


We  now  have  a set  of  simultaneous  linear  equations  given  by 


N 


[3.30] 


Since  the  matrix  w is  typically  very  large  and  sparse,  the  equations  are  normally  highly 
underdetermined  (M«N),  and  errors  in  the  data  cause  them  to  be  inconsistent,  the 
solution  is  best  done  iteratively.  Since  the  calculation  of  wy  as  given  in  [3.29]  adds  a 
significant  amount  of  computational  complexity  to  the  problem  it  is  typically 
approximated  as  follows.  If  the  centroid  of  element  is  contained  in  the  projection 
path  ^ , then  Wij  is  given  a value  of  one,  otherwise  it  is  set  tozero. 

The  ART  reconstruction  process  is  carried  out  as  follows.  Let  Ui‘’  represent  the  q‘*’ 
estimate  of  nj  and  represent  the  q'*’  estimated  value  of  the  projection  data  point  Pj. 

1)  An  initial  estimate  is  given  to  n.  Typically  nj*^=  iio  or  ni'^^O  for  all  i is  used. 

2)  Initialize  the  projection  angle  to  the  first  projection  angle. 

3)  Initialize  the  projection  data  point  to  the  first  projection  data  point.  Pi. 

4)  Compute  the  difference  between  the  actual  projection  data  and  the  estimate  from 
the  sum  of  the  reconstruction  elements  as  given  in  [3.29]. 

5)  This  difference  is  then  divided  among  the  reconstruction  elements  which 
contributed  to  it  in  [3.29],  in  proportion  to  their  weight,  Wy,  by  adding  it  to  the 
previous  values. 


6)  Move  to  the  next  projection  data  point  and  repeat  from  step  4. 
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7)  After  all  projection  data  points  for  the  present  projection  angle  have  been  used, 
move  to  the  next  projection  angle  and  repeat  from  step  3. 

The  seven  steps  described  above  represent  the  first  iteration  of  the  reconstruction  process 
and  result  in  Uj’ . The  next  iteration  is  then  started  over  from  step  2.  Note  that  the 
measured  projection  data  is  typically  in  the  form  of  phase  shift  rather  than  absolute  phase, 
so  that  we  are  really  reconstructing  an  index  difference  , - n{r))  rather  than  n(r) . 

However  , the  above  description  still  applies  with  no  loss  of  generality  if  (n^  - n(r))  is 
substituted  for  n(r) . 

To  determine  when  the  algorithm  has  converged  to  an  acceptable  solution,  it  is 
necessary  to  examine  some  measure  of  the  discrepancy  between  the  measured  and 
calculated  projection  data.  Gordon  et  al.  (1970)  and  Gordon  (1974)  have  devised  a 
number  of  such  measures. 

ART  Variations 

There  are  a number  of  variations  on  ART  which  may  give  differing  results 
depending  upon  the  application.  Multiplicative  ART  or  MART  follows  the  same  basic 
procedure  as  ART  with  the  exception  of  the  iterative  error  update.  In  MART  the  error 
between  the  calculated  and  measured  projection  data  is  applied  to  the  reconstruction 
elements  as  a multiplicative  factor  as  follows. 


nr'  = 


yPh 


[3.31] 


Another  variation  of  the  technique  is  to  set  any  reconstructed  index  value  equal  to  zero 
when  the  above  calculations  result  in  a negative  value.  This  is  a form  of  constrained 
ART.  Other  constraints,  which  may  be  unique  to  a specific  application,  could  also  be 
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applied  and  this  is  one  of  the  important  features  of  the  iterative  methods  which  is 
exploited  when  the  projection  data  is  limited  by  an  obscuration.  Cha  and  Sun 
(1989,1990)  have  addressed  this  problem  in  significant  detail. 

Tomographic  Reconstruction  Noise  Analysis 
Convolutional  Back-Projection  Closed  Form  Solution 

To  gain  insight  into  the  process  by  which  noise  contained  in  the  projection  data  is 
carried  through  the  tomographic  reconstruction  process,  we  will  examine  the 
convolutional  back-projection  algorithm  for  which  an  analytical  solution  exists. 
(Kak&Slaney,  1988)  We  will  also  consider  other  reconstruction  algorithms,  such  as 
iterative  methods,  for  which  an  analytical  solution  carmot  be  obtained.  In  these  cases 
there  is  empirical  data  available  in  the  literature  which  has  been  obtained  from  monte- 
carlo  simulations  and  measurements.  In  order  to  develop  an  analytical  solution  for  the 
convolutional  back-projection  algorithm,  we  must  first  make  some  fundamental 
assumptions  regarding  the  noise  on  the  projection  data.  For  the  following  analysis,  we 
will  assume  additive,  wide  sense  stationary,  zero-mean,  white  noise.  Under  these 

assumptions,  the  noisy  projection  data,  Pg{x)  ,can  be  written  in  terms  of  the  noise  free 
projection,  Pg{x) , and  a noise  term,  rig{x)  . 

Pe{x)  = Pg{x)  + ng{x)  [3.32] 

In  theory,  the  white  noise  assumption  on  the  projected  phase  shift  measurements 
allows  for  negative  values  of  integrated  phase  shift,  which  would  only  be  possible  for 
fluid  densities  greater  than  the  ambient  density.  Although  this  is  not  physically  possible 
in  a pure  aerodynamic  flow  measurement,  no  attempt  will  be  made  to  nullify  negative 
values,  which  can  occur  in  the  measured  phase  shifts  due  to  thermal  noise  on  the 
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wavefront  sensor  detector  pixels.  It  may  be  possible  to  improve  the  accuracy  of  the 
tomographic  reconstruction  by  performing  this  nullification  of  the  projection  data, 
although  the  noise  analysis  would  become  intractable. 

Under  the  assumptions  that  iig{x)  is  a stationary  random  process  which  is 
uncorrelated  from  one  sample  to  the  next  as  well  as  from  one  projection  angle  to  the  next 
we  can  write, 


If  Sg{^  is  the  Fourier  transform  of  the  noisy  projection  data,  Pg{x) , then  the  filtered 
projection  data  can  be  written  as 


where  G(^  is  the  smoothing  filter  including  the  finite  detector  width  sampling  effects. 
Back-projection  of  the  filtered  data  yields. 


Recall  that  the  projection  data  contains  a deterministic  part , Sg{(^ , as  well  as  a random 
part,  Ng{(^  , due  to  the  projection  data  noise  where. 


[3.33] 


Q;(x)=  , 


[3.34] 


[3.35] 


0 


or,  in  terms  of  the  projection  data  in  the  Fourier  domain. 


[3.36] 


00 


-00 


[3.37] 
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We  now  wish  to  find  the  statistics  of  the  reconstructed  data,  / (;r,  . The  mean  can  be 

found  as  follows. 


E 


f(x,y) 


0 -00 


[3.38] 


The  only  random  part  is  Ng{^  . 


£[N,(f)]  = E 


00  00 

j ng(x)  = j E[ng(x)]  = 0 


[3.39] 


So  that, 


f{x,y) 


= I J 5,(«l«|G(^)< 


j‘2/t4(xcos0+ysin 


^^d^dO 


[3.40] 


^ -*  0 -00 

which  is  just  the  noise  free  reconstruction.  The  variance  of  the  reconstructed  data  can  be 
written  as 


<^\{x,y)  = E 


f{x,y)-E 


f(x,y) 


2 


[3.41] 


Since  the  deterministic  term  can  be  removed  from  the  expectation  operation  for  both 
terms,  it  cancels,  and  since  £'[N^((^)]  = 0 this  leaves. 


al(x,y)  = E 


= E 


0 —00 

!|  j A'.,  (^,  )|« , 1G(« , dS,  I X 


0 —00 


[3.42] 


The  expectation  can  be  moved  inside  the  integrand  and  since  Ng{^  is  the  only  random 


part. 
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= j j j {QK  (4)]  1^,1141  d^,d^,  de,  de,  ■ [3.43] 

0 0 -00-00 

We  can  rewrite  E^N  N*g^{^2)\  as  follows: 


NeS^x)N*eS^2)]=  I J^K(^i)  • [3.44] 


Under  the  assumption  of  uncorrelated  noise,  equation  [3.33]  holds  so  that, 


^'.,(^,)'v;(^2)]=  -x,)S(g,  -e,)e’^’^^'’'-^'’'Ux,dx,  .[345] 


By  the  sifting  property  of  the  delta  function 


E[N.,(i,)  Nl  (f, )]  = I iV.  S(e,  - e, ) 

—00 

00 

= J(6>,  -6»2)  dx,  . 


[3.46] 


00 

Now  since  J j can  write, 

—00 


E[N.,{i,)Nl^Xi,)]  = N„S(e,-e,)S{^,-i,)  . [3.47] 

So  that  we  have, 

0 0 -00-00 

Using  the  sifting  property  again,  we  are  left  with 

a\{x,y)  = Nj\]%\^\G{^,)fd^,de,=7vN,]\^^\G{^fd^  • [3-49] 

0 -00  -QO 

Since  the  original  projection  data  was  zero  mean,  its  variance  is  E^nl  (x)j  - N^.  So  that 
the  noise  power  is  increased  by  a factor  K where. 
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[3.50] 


Recall  that  G{^  contains  rect  -^j  = rect[^a)  always  and  may  also  include  additional 

filtering,  such  as  finite  detector  width  effects.  Also  note  that  in  the  discrete 
implementation,  the  integration  of  the  continuous  function  over  theta  in  [3.49]  becomes  a 
summation  of  M discrete  projections  and  we  must  divide  the  reconstruction  by  kg  to  get 
back  to  units  of  refractive  index.  The  discrete  version  of  [3.49]  is  then  as  follows. 


^^0 

klM 


M “ 
m=  \ —00 


[3.51] 


This  relationship  can  also  be  written  in  terms  of  a factor  K,  which  signifies  the  increase  in 
the  projection  variance  manifest  in  the  reconstruction,  and  can  be  written  as  follows: 


o -w 

Emperical  Data  for  Iterative  Algorithms 

As  discussed  above,  when  the  projection  data  is  limited  by  opaque  objects  within 
the  index  field,  or  by  projection  angle  limitations,  iterative  reconstruction  algorithms 
must  be  used.  The  effect  of  projection  data  noise  on  reconstruction  accuracy  using 
iterative  algorithms  has  been  studied  by  a number  of  authors.  Bahl  and  Liburdy  (1991), 
present  the  results  of  reconstruction  using  the  ART  iterative  algorithm  in  comparison  to 
convolutional  back-projection  as  well  as  several  others.  Verhoeven  (1993),  presents 
similar  work  using  five  different  iterative  algorithms  which  are  adapted  versions  of  ART 
and  MART.  In  these  studies,  an  analytical  index  profile  is  assumed  and  simulated  noise 
of  differing  levels  is  added  to  the  numerically  computed  projection  data.  The  index 
profile  is  then  reconstructed  using  the  algorithm  of  interest  and  reconstruction  errors  are 
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computed.  In  Verhoeven's  work,  average  reconstruction  errors  of  less  than  5%  are 
achieved,  using  the  ART  algorithm  and  19  projection  angles  spaced  over  90  degrees,  with 
a projection  noise  levels  of  20  dB  SNR  or  better.  When  the  SNR  was  increased  to  30  dB, 
the  average  reconstruction  error  was  nearly  equal  to  the  error  obtained  with  no  noise. 

Some  general  conclusions  can  be  drawn  from  this  work.  An  important  result  is 
that  the  reconstruction  accuracy  is  dependent  upon  the  characteristics  of  the  index  profile 
and  upon  the  number  of  iterations.  When  the  projection  data  contains  noise,  the 
minimum  reconstruction  error  is  obtained  with  fewer  iterations  than  when  noise  is  not 
present.  There  is  also  significant  variation  in  reconstruction  accuracy  between  algorithms 
depending  upon  the  amount  of  projection  data  that  is  available. 

For  a specific  application  of  one  of  the  iterative  reconstruction  algorithms  it  may 
be  prudent  to  perform  a similar  analysis  under  the  specific  conditions  of  that  application. 
This  type  of  analysis  would  yield  a more  accurate  determination  of  the  expected 
reconstruction  accuracy  as  a function  of  projection  noise  level.  This  will  be  done  in  the 
next  section 

Design  Example 

Although  the  tomographic  reconstruction  process  is  not  the  focus  of  this  work,  it  is 
important  to  understand  its  concepts  and  the  requirements  imposed  by  it  upon  the 
projection  data.  To  this  end,  this  chapter  has  provided  the  background  necessary  to 
understand  the  guidelines  used  to  establish  requirements  on  the  projection  data.  Using 
the  supersonic  conical  flow  model  as  a design  example,  we  can  now  apply  some  of  the 
results  of  this  chapter.  The  key  parameters  of  interest  in  this  example  are  as  follows. 


1)  Test  region  Size:  >50mm 


65 


2)  Spatial  Resolution:  < \mm 

3)  Wavelength:  532nm 
Projection  Data  Sampling 

The  projection  data  sample  spacing  can  be  easily  determined  by  applying  the 
Nyquist  theorem.  Since  a spatial  resolution  of  no  more  than  \mm  is  required,  the  sample 
spacing  must  be  no  larger  than  0.5mm,  which  is  half  the  spatial  frequency  resolution 
period.  However,  if  we  consider  the  low  pass  sampling  effects  of  the  large  sampling 
aperture  of  a 100%  duty  factor  lens  array,  we  must  apply  the  filter  function  described  by 
[3.26].  At  the  3dB  point,  the  argument  of  the  sine  function  takes  on  a value  of  1.3823 
which  can  be  applied  to  determine  the  sample  spacing  as  follows: 


This  results  in  a minimum  of  1 13  samples  for  each  slice  of  the  50mm  test  region,  which 
is  a 1 13  by  1 13  array  of  samples  of  a 1 13X1 13  Hartmann  lens  array. 

Required  Number  of  Projection  Angles 

Applying  [3.9]  to  the  design  example  is  straightforward. 


As  stated  previously,  the  relationship  given  by  [3.9]  is  intended  as  a very  rough  rule  of 
thumb.  It  is  based  upon  the  spatial  resolution  requirement  at  the  edges  of  the  test  region 
and  results  in  over-sampling  toward  the  center  of  the  test  region.  Obviously  a 
requirement  for  179  projection  angles  spread  over  180  degrees  would  result  in  a 
prohibitive  system  complexity  for  most  high-speed  flow  measurement  applications.  As  a 
result,  it  would  now  be  important  to  examine  the  requirements  of  an  iterative 


= 1.3823^  a = 


2 


2(1.3823)  _ 2(1.3823)  _ 1.3823  _ 1.3823 

^3dB  3dB  ^ fidB  ^ \mm 


= O.AAmm . [3.53] 


Ttd  jz{50mm) 

= — V = 178.5 

2 a 2{0.AAmmj 


[3.54] 


reconstruction  technique  by  using  simulated  projection  data  and  examining  the 
reconstruction  performance  as  the  number  of  projection  angles  is  varied. 
Reconstruction  Noise  Variance 
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The  relationships  given  in  [3.51]  and  [3.52]  can  now  be  applied  to  give  some 
insight  into  how  noise  on  the  projection  data  manifests  itself  in  the  reconstruction.  For  a 
100%  duty  factor  lens  array,  the  sampling  filter  given  by  [3.26]  becomes. 


sin(;r^a) 

K^a 


[3.55] 


Equation  [3.37]  is  now. 


Mia 


sm{7i^a) 

n^a 


[3.56] 


Substituting  for  M=18  projection  angles,  A,=532nm,  a=0.44  mm  and  integrating 

Q 

numerically  yields  a value  for  K of  approximately  1.339  x 10'  . Assuming  a projection 
data  standard  deviation  of  10  degrees  of  phase,  would  result  in  a reconstructed  index  of 
refraction  noise  variance  as  follows: 


cr^  =yjKa^p  = 1.339x  10' 


10—1  =2.02x10'^ 
1,  360  j 


[3.57] 


Note  from  Figure  2-2  that,  for  the  given  design  example,  this  corresponds  to  a 
rather  large  RMS  error  of  18.4%  of  the  peak  index  difference. 

Tomographic  Reconstruction 

For  the  sake  of  completeness,  a reconstruction  of  the  index  profile  of  the  design 


example  was  made  using  both  the  convolutional  back-projection  (CBP)  method  and  the 
iterative  ART  method.  To  allow  the  use  of  CBP,  a constant  index  value  equal  to  the 
boundary  value  was  assumed  for  the  region  within  the  cone  body.  This  allowed 
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projection  data  to  be  calculated  and  reconstructions  to  be  made  without  regard  to  the  data 
obscuration  problem.  The  reconstructions  were  made  using  projection  data  calculated  for 
18  projection  angles  spaced  by  10  degrees,  although  the  axial  symmetry  of  this  example 
resulted  in  the  same  data  for  every  angle.  Figure  3-6  shows  the  reconstruction  from  the 
CBP  method  and  Figure  3-7  shows  a radial  cut  of  this  data  plotted  along  with  the  original 
index  profile.  Figure  3-8  shows  the  reconstruction  from  the  iterative  ART  method  and 
Figure  3-9  shows  a radial  cut  of  this  data  plotted  along  with  the  original  index  profile. 
Since  the  index  profile  was  known  a priori,  it  was  possible  to  calculate  the  reconstruction 
errors,  and  the  ART  reconstruction  was  terminated  when  the  RMS  error  began  to 
increase.  The  CBP  reconstruction  had  a RMS  error  of  1 .03x10'^  and  a peak  error  of 
4.65x10'^.  The  ART  method  began  to  diverge  after  the  7th  iteration  at  which  point  the 
RMS  error  was  2.50x10'^  and  the  peak  error  was  2.65x10'^.  These  RMS  errors 
correspond  to  9.5%  and  2.3%  of  the  peak  index  difference  for  the  CBP  and  ART 
reconstructions,  respectively.  Note  that  the  index  profile  including  the  sharp  edge  at  the 
shock  boundary  was  reconstructed  accurately  with  much  less  than  the  number  of 
projections  called  for  by  [3.9]. 

In  order  to  verify  the  noise  analysis  and  the  results  of  equation  [3.51], 
reconstructions  were  performed  on  the  design  example  index  profile  with  added 
projection  data  noise.  Zero-mean,  Gaussian  noise  was  added  to  the  projection  data, 
which  was  uncorrelated  from  sample  to  sample  as  well  as  from  one  projection  angle  to 
another,  and  the  iterative  ART  reconstruction  was  performed  as  well  as  the  convolutional 
back-projection  algorithm.  Figure  3-10  shows  how  the  RMS  error  increased  as  a function 
of  the  standard  deviation  of  the  projection  data  noise,  as  a percentage  of  the  peak  index 
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Figure  3-6:  CBP  reconstruction  of  design  example  index  profile  difference 


Figure  3-7:  CBP  reconstruction,  comparison  to  original  index  difference 


Figure  3-8:  ART  reconstruction  of  design  example  index  profile  difference 


(Do-n)  X 10-'' 


Figure  3-9:  ART  reconstruction,  comparison  to  original  index  difference 
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difference  (1.0799x10''^).  The  reconstruction  standard  deviation  of  the  ART 
reconstruction  is  well  below  that  resulting  from  the  CBP  algorithm,  and  both  are  better 
than  the  error  predicted  by  equation  [3.51].  Projection  data  noise  standard  deviations 
from  5 to  60  degrees  were  used  for  a projection  slice  over  which  the  maximum  phase 
shift  was  2715  degrees.  These  values  of  noise  standard  deviation  correspond  to  0.2  to  2.2 
percent  of  the  maximum  phase  shift  or  27.3  to  16.6  dB  in  terms  of  signal  to  noise  ratio. 

As  seen  in  Figure  3-10,  the  reconstructions  contain  errors  with  no  projection  data  noise, 
which  are  due  to  the  reconstruction  itself 


Figure  3-10;  RMS  reconstruction  errors 
(ART-solid  line,  CBP-dashed  line,  calculated  by  [3.51]-dash/dotted  line) 
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Conclusions 

The  intent  of  this  chapter  has  been  to  provide  enough  background  on  the 
tomographic  reconstruction  process  so  that  the  discrete  sampling  and  noise  requirements 
for  the  projection  data  could  be  established.  One  important  result,  which  should  be 
stressed,  is  that  the  reconstruction  accuracy  and  projection  data  requirements  are  highly 
dependent  upon  the  index  profile  being  measured  and  on  the  details  of  the  reconstruction 
algorithm.  The  reconstruction  algorithm  must  be  tailored  to  the  expected  index  profile 
and  to  the  measurement  system  specifics  such  as  obscurations.  The  results  of  the  design 
example  analysis  above  indicate  that  a good  reconstruction  of  this  tjqie  of  index  profile  is 
possible  under  the  following  conditions. 

1)  Projection  Data  Sample  Spacing:  >113  points  spaced  at  <0.44  mm. 

2)  Number  of  Projection  Angles:  1 8,  evenly  spaced  over  1 80  degrees. 

3)  Projection  Data  Noise:  Standard  Deviation  < 30°  or 

SNR>19.6  dB 

The  18  projection  angles  used  here  are  well  below  the  number  called  for  by 
equation  [3.54].  However,  18  projection  angles  will  still  result  in  considerable  cost  and 
complexity  in  terms  of  the  optical  system.  Building  fewer  projection  angles  could  still 
produce  useful  results  in  some  measurement  applications,  although  it  would  be  highly 
dependent  upon  the  measurement  being  performed.  Measurements  which  have  nearly 
axial  symmetry  and  for  which  the  spatial  resolution  requirement  can  be  relaxed  near  the 
boundaries  of  the  test  region  would  lend  themselves  to  a system  with  fewer  projections. 
Likewise,  asymmetric  test  regions  for  which  the  spatial  resolution  is  critical  throughout 
the  test  region  would  require  more  than  the  1 8 projections  used  here.  For  this  reason,  it  is 
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important  to  test  the  reconstruction  algorithm  using  simulated  projection  data,  which  is 
representative  of  the  type  of  measurement  being  performed. 


CHAPTER  4 

WAVEFRONT  PHASE  RECONSTRUCTION 
Problem  Statement 

The  tomographic  reconstruction  process  uses  integrated  phase  shift  data  to  reconstruct 
the  fluid  density  in  the  test  region.  Since  the  wavefront  sensor  measures  the  phase 
gradient,  an  intermediate  stage  is  required,  in  which  the  phase  shift  is  reconstructed  from 
the  phase  gradient  data.  Recall  that  the  measured  spot  deviations  represent  the  localized 
wavefront  gradient  averaged  across  the  sub-aperture  associated  with  each  spot.  The 
wavefront  reconstruction  problem  becomes  one  of  fitting  this  data  by  means  of  a least 
squares  or  some  other  criteria  to  an  estimate  of  wavefront  phase. 

Background 

The  problem  of  reconstructing  the  wavefront  phase  shift  from  a set  of  wavefront 
sensor  measurements  has  been  studied  extensively  for  the  case  of  atmospheric  distortion 
and  optimum  solutions  based  on  various  techniques  have  been  derived.  The  techniques 
used  fall  into  one  of  two  classes  consisting  of  zonal  or  modal  reconstructions. 

Zonal  methods 

The  first  class  consists  of  zonal  methods  in  which  the  estimate  is  iteratively  matched 
to  the  measured  data  for  each  spatial  zone.  Zonal  techniques  were  used  by  Hudgin 
(1977)  as  well  as  by  Fried  (1977)  in  deriving  an  optimal  estimate  of  the  wavefront  phase 
from  a set  of  wavefront  sensor  measurements  defined  over  a square  array  and  using  a 
least  squares  solution.  Both  authors  showed  that  the  mean  square  wavefront  error  due  to 
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measurement  noise  depends  on  the  number  of  measurement  points,  N,  according  to  A + 
ln(N).  Noll  (1978)  derived  the  same  result  analytically.  Hunt  (1979)  formulated  the  same 


problem  using  matrix  methods  and  developed  a more  efficient  reconstruction  scheme  in 
terms  of  its  rate  of  convergence. 

Modal  Methods 

Cubalchini  (1979)  examined  the  reconstruction  problem  from  a modal  point  of  view, 
which  is  the  second  class  of  techniques  used.  In  this  method  the  estimate  is  based  on  a 
set  of  orthogonal  basis  functions,  which  can  be  represented  by  a set  of  coefficients.  As 
with  most  of  the  research  conducted,  Zemike  polynomials  were  used  as  the  set  of  basis 
functions  because  they  are  well  matched  to  atmospheric  aberrations  and  astronomical 
telescope  applications.  Southwell  (1980)  and  Herrmann  (1980)  examined  the 
reconstruction  problem  from  both  points  of  view  and  made  comparisons  between  them. 
Southwell  (1980)  showed  that  zonal  methods  are  superior  in  terms  of  noise  propagation 
as  well  as  being  more  computationally  efficient  and  that  the  superiority  is  enhanced  with 
an  increasing  number  of  measurements. 

Choice  of  Method 

Although  prior  research  has  shown  that  the  modal  techniques  are  superior  to  zonal 
ones  for  the  atmospheric  distortion  correction  application,  this  is  only  true  in  cases  where 
a suitable  set  of  orthogonal  basis  functions  is  available.  The  basis  functions  must  be  well 
matched  to  the  expected  reconstructions,  as  the  Zemike  polynomials  were  in  this 
application.  In  order  to  make  the  reconstmctions  more  useful  for  a broader  class  of 
measurement  applications,  and  since  computational  efficiency  is  of  less  concern  for  this 
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non-real-time  application,  we  have  chosen  to  use  the  more  general  zonal  reconstruction 
approach. 

In  our  application,  the  wavefront  sensor  yields  N X N wavefront  slope 
measurements  in  one  direction  and  another  N X N wavefront  slope  measurements  in  the 
orthogonal  direction.  If  the  wavefront  sensor  lies  in  the  x-y  plane  and  the  wavefront  tilt 
angle  in  the  x-z  plane  for  lenslet  i,j  is  Qk.  . then  the  focussed  spot  x directed  deviation 
5x.  j is  simply  related  by  the  lenslet  focal  length  / , as 

dx^  j = / tan  6x^  j s / 6x.  j . [4.1] 


The  wavefront  tilt  angle  &x.  . can  be  related  to  the  phase  slope  (jt'x^  j across  a 


lenslet  of  aperture  diameter  a,  as  follows: 
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, , In  _ 2 ;r  ^ ( rad 

(p  x- = tan  6fc,. s Ox, , 

^ '■'  A '''  A '•'  yum. 


[4.2] 


We  can  represent  the  2N^  phase  slope  measurements  using  a grid  of  (N+1)^  phase  data 
points  as  in  Figure  4-2.  The  phase  slope  can  be  written  as  a first  order  fit  to  the  phase 

9 9 

grid  points,  generating  2N  equations  with  (N+1)  unknowns  as  follows: 


[4.3] 


Since  we  are  only  interested  in  the  relative  phase  shift  over  the  reconstruction  grid,  we 
can  use  an  arbitrary  boundary  condition  for  the  absolute  phase  such  as  (px,  j = (py,  ^ = 0 . 
By  writing  the  2N^  equations  in  matrix  form  we  can  efficiently  find  the  least  squares 
solution  using  singular  value  decomposition  methods.  The  equations  for  the  x-directed 
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phase  slopes  and  the  y-directed  phase  slopes  can  be  put  into  a similar  form  and  appended. 
This  will  yield  one  matrix  equation  in  which  a sparse  [2N^  x (N+1)  matrix  multiplied 
by  a [(N+1)  ^ x 1]  vector  of  unknown  phase  values  is  equated  to  a [2N^  x 1]  vector  of 
measured  x and  y directed  phase  slopes. 
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The  least  squares  solution  to  this  equation  will  be  a first  order  fit  of  the  relative  phase 
shift  to  the  measured  phase  gradient  data.  Since  the  wavefront  sensor  measures  the  phase 
slope  only,  any  information  beyond  the  first  order  is  lost  in  the  measurement  and  no 


78 


benefit  is  gained  by  using  a higher  order  fit.  As  long  as  the  phase  slope  data  is  sampled  at 
or  above  the  Nyquist  criteria,  this  will  be  sufficient  for  a good  reconstruction.  A 
wavefront  phase  reconstruction  algorithm  of  this  type  has  been  developed  in  support  of 
this  research  (Gunia,  1999).  This  algorithm  was  developed  using  Matlab,  and  performs 
the  least  squares  fit  as  outlined  above  using  singular  value  decomposition.  The 
reconstruction  is  done  in  a piecewise  fashion  in  order  to  speed  up  the  reconstruction 
process  to  reasonable  time  periods.  One  important  result  of  Gunia's  work  was  to 
demonstrate  the  importance  of  the  joining  process  used  in  the  piecewise  reconstruction. 

By  using  the  average  of  the  boundary  values  in  joining  adjacent  pieces  of  the 
reconstruction,  the  overall  error  can  be  reduced.  This  is  the  approach  used  here. 

Reconstruction  Algorithm  Performance 
In  order  to  examine  the  expected  accuracy  of  the  least  squares  reconstruction 
algorithm,  simulated  noisy  wavefront  sensor  data  was  used  with  various  overall  and  inter- 
algorithm reconstruction  sizes.  Gaussian  noise  with  standard  deviations  ranging  from  10 
to  30  dB  below  the  peak  spot  deviation  was  added  to  the  ideal  spot  deviation  data  from  a 
linearly  tilted  wavefront.  The  reconstruction  was  then  performed  and  the  errors  were 
calculated.  Using  the  errors  from  all  points  in  the  reconstruction  grid,  the  standard 
deviation,  or  RMS  error,  of  the  reconstruction  was  calculated.  A number  of  simulations 
were  performed  to  determine  that  the  reconstruction  error  was  independent  of  the 
wavefront  tilt  value.  Ten  simulations  for  each  noise  level  were  performed  and  the 
average  standard  deviation  was  determined.  An  overall  reconstruction  sizes  ofll4X  114 
was  used  with  inter- algorithm  reconstruction  grids  of  3X3,  5X5,  7X7,  9X9,  1 1X1 1, 
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13X13,  and  15X15.  Figure  4-3  shows  the  results  for  the  1 14X1 14  reconstruction.  Note 
from  the  figure  that  as  the  inter-algorithm  grid  size  is  increased,  the  reconstruction  error 


1 14  X 114  Reconstruction  Noise  Standard  Deviation 
Block  Size  (3-15) 


Figure  4-3:  114X114  reconstruction  accuracy 


is  decreased.  When  the  noise  level  is  20  dB  below  the  peak  deflection  limit,  a 
reconstruction  grid  size  of  3 results  in  a reconstruction  RMS  error  of  28.3  degrees  while  a 
grid  size  of  15  results  in  15  degrees  of  RMS  error.  As  the  noise  level  approaches  the  -30 
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dB  level,  there  is  little  difference  between  the  various  grid  sizes.  However,  the  use  of 
larger  piecewise  reconstruction  grid  sizes  becomes  more  important  as  the  noise  level  is 
increased  toward  -10  dB.  Also  note  that  the  reconstruction  errors  decrease  as  the  overall 
size  of  the  reconstruction  is  decreased.  This  is  in  agreement  with  the  research  of  Hudgin 
(1977)  , Fried  (1977)  , and  Noll  (1978). 

Design  Example 

Using  the  supersonic  conical  flow  design  example,  a wavefront  reconstruction 
was  done  on  simulated  noise-free  Hartmann  sensor  data.  This  wavefront  sensor  data  has 
discontinuities  due  to  the  shock  and  the  obscuration  by  the  cone.  The  simulated  phase 
shift  data  is  shown  in  Figure  4-4.  The  wavefront  phase  reconstruction  program 


Simulated  Wavefront  Phase  Shift  Data 


Figure  4-4:  Simulated  wavefront  phase  shift  data 
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developed  by  Gunia  (1999)  was  not  intended  for  data  with  discontinuities  or 
obscurations.  The  wavefront  sensor  data  for  the  obscured  region  was  zeroed  in  order  to 
allow  the  reconstruction  to  proceed  with  a full  array  of  data.  As  a result  of  the  resulting 
discontinuities,  the  wavefront  phase  reconstruction  contained  significant  errors.  More 
accurate  wavefront  reconstructions  are  possible  by  using  more  complex  reconstruction 
algorithms  tailored  to  discontinuous  data.  It  is  likely  that  the  principles  developed  by 
Cha  and  Sun  (1990)  or  Vest  and  Prikryl  (1984)  for  tomographic  reconstructions  with 
obscurations  and  discontinuities  would  be  directly  applicable  to  the  wavefront 
reconstruction  problem.  Both  of  these  authors  used  available  a priori  information  and 
constraints  along  with  iterative  solutions  to  subvert  the  problems  encountered  with 
discontinuous  and  obscured  data.  Since  the  development  of  such  an  algorithm  was 
outside  the  scope  of  this  work,  the  wavefront  reconstruction  program,  developed  for 
continuous  data,  was  used  with  a small  reconstruction  grid  size.  The  reconstruction  was 
performed  for  the  sake  of  completeness,  and  although  the  resulting  errors  were  large,  the 
reconstruction  takes  on  the  basic  form  of  the  simulated  data.  Figure  4-5  shows  the 
reconstructed  wavefront  phase  using  a reconstruction  grid  size  of  4.  The  RMS  error  was 
667  degrees,  corresponding  to  a peak  SNR  of  around  9 dB.  This  is  obviously  well  below 
the  acceptable  level  of  data  quality  to  allow  an  accurate  tomographic  reconstruction. 

Summary  and  Conclusions 

In  summary,  a zonal  minimum-least-squares  fit  algorithm,  as  detailed  above,  will 
be  used  to  reconstruct  the  wavefront  phase  from  the  phase  gradient  data.  The  parameters 
of  this  algorithm  must  be  set  such  that  the  resulting  wavefront  phase  reconstruction  errors 
are  less  than  the  maximum  allowable  error  required  by  the  tomographic  reconstruction 


algorithm.  The  results  of  the  previous  chapter  indicated  that  for  the  supersonic  conical 
flow  design  example,  accurate  tomographic  reconstructions  are  made  with  1 -sigma 
projected  phase  errors  on  the  order  of  30  degrees.  The  wavefront  reconstruction  on  the 
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Reconstructed  Wavefront  Phase  Shift  Data 


6000 


5000.  ’ 


Figure  4-5  : Wavefront  reconstruction  (4X4  reconstruction  grid  size) 


simulated  design  example  resulted  in  errors  on  the  order  of  700  degrees  due  to  the  nature 
of  the  algorithm  alone.  For  more  realistic  data,  which  may  contain  shock  discontinuities 
and  opaque  obscurations,  a more  robust  wavefront  phase  reconstruction  program  is 
required.  This  is  an  area  where  more  research  is  needed. 


CHAPTER  5 

SPOT  DEFLECTION  MEASUREMENT 


Problem  Statement 

As  described  in  the  previous  chapter,  in  order  to  generate  phase  data  with  the 
precision  required  by  the  tomographic  reconstruction  algorithm,  the  spot  deflection 
measurement  must  have  a dynamic  range  of  at  least  20  dB.  Specifically,  the  RMS  noise 
level  of  the  spot  deflection  measurement  must  he  less  than  or  equal  to  1/1 of  the 
maximum  anticipated  deflection.  The  characteristics  of  the  spot  deflection  measurements 
are  determined  by  the  lens  array  design,  the  detector  array  geometry,  and  the  spot 
position  measurement  strategy.  For  the  purposes  of  this  work  we  will  assume  that  we 
have  a square  grid  of  lenslets  in  front  of  a square  grid  of  detector  pixels.  In  order  to 
maintain  generality,  we  will  allow  for  a stage  of  magnification,  between  the  lens  array 
and  detector  array,  for  matching  purposes.  The  projected  footprint  of  each  lenslet  will 
subtend  a group  of  N'  by  N'  detector  pixels  of  which  a subset  of  N by  N pixels,  centered 
on  the  focussed  spot,  are  used  to  determine  the  position  of  the  spot.  Centroiding  or 
curve-fit  methods  can  be  used  to  calculate  the  spot  position  to  sub-pixel  accuracy.  Figure 
5-1  shows  an  example  of  this  spot  measurement  geometry.  Assuming  that  the  size  of 
each  detector  sub-array,  N'  is  set  such  that  the  maximum  anticipated  deflection  will  fall  at 
the  edge  of  the  sub-array,  then  our  20  dB  dynamic  range  requirement  requires  that  the 
RMS  error  be  less  than  1/100^’’  of  this  amount  of  maximum  deflection.  We  will  assume 
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that  the  focussed  spot  intensity  profile  is  diffraction  limited.  The  two  types  of  errors 
which  are  the  primary  contributors  to  the  overall  wavefront  tilt  measurement  accuracy  are 
errors  due  to  the  discrete  nature  of  the  centroiding  or  curve  fit  algorithm  and  errors  due 


to  the  detector  measurement  noise.  A tradeoff  must  be  made  to  determine  how  many 
detector  pixels  are  used  in  the  centroid  or  curve  fit  calculations.  Too  many  pixels  used 
will  result  in  increased  noise  errors  while  too  few  pixels  will  cause  large  algorithm  errors. 
Secondary  sources  of  error  are  those  due  to  lenslet  focal  length  uncertainty  and 
misalignment  between  the  lens  array  and  detector  array. 


Background 

There  is  a broad  base  of  information  available  concerning  the  accuracy  of  centroid 
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measurement  using  a CCD  array,  although  most  of  the  research  was  not  done  with  a 
Hartmann  sensor  in  mind.  Tyler  and  Fried  (1982)  developed  the  image  position  error 
associated  with  a quadrant  detector  as  a function  of  the  signal  to  noise  voltage  ratio. 
Winick  (1986)  determined  the  Cramer-Rao  lower  bounds  for  this  measurement  using  one 
or  two  dimensional  detector  arrays  and  Down  (1992)  has  analyzed  the  image  position 
error  for  focal  plane  arrays  by  assuming  poisson  statistics  for  object  and  background  shot 
noise  and  detector  dark  noise.  Grossman  and  Emmons  (1984)  examined  the  problem  of 
optimizing  a detector  array  for  point  source  tracking  using  odd-point  centroids  and 
quadratic  curve  fits.  Much  of  this  work  can  be  applied  to  the  Shack-Hartmann  sensor 
design. 

Spot  Position  Measurement  Errors 

Algorithm  Error 

The  problem  of  finding  the  centroid  or  curve  fitting  to  find  the  spot  position  uses 
discretely  sampled  data  in  the  form  of  integrated  intensities  of  discrete  detector  array 
pixels  to  estimate  the  continuous  function  of  spot  position.  A large  number  of  samples 
(pixels)  are  required  to  accurately  estimate  the  position  across  the  full  range  of  the 
measurement.  As  a result,  in  most  practical  applications  there  will  be  significant 
algorithm  errors  caused  by  the  use  of  a limited  amount  of  data  in  the  centroid  or  curve  fit 
calculations.  Fortunately,  these  errors  are  deterministic  and  they  can  be  calibrated  out  to 
the  extent  that  they  can  be  accurately  modeled. 
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In  a typical  spot  centroid  measurement  strategy,  such  as  that  shown  in  Figure  5-1,  a 
coarse  position  estimate  is  required  to  find  the  pixel  with  the  peak  output.  Next,  a subset 
of  N-by-N  pixels  is  used  to  find  the  centroid,  with  typical  values  for  N being  in  the  range 
of  2 to  7.  The  centroid  is  found  in  one  direction  by  summing  the  N rows  (or  columns)  in 
the  orthogonal  direction,  and  calculating  an  N-point  geometric  centroid  of  the  detector 
intensity  values.  The  summation  is  then  reversed  to  find  the  N point  centroid  in  the  other 
direction.  An  N odd-point  centroid,  would  be  calculated  with  respect  to  the  center  of  the 
peak  valued  pixel p,q  as  follows.  Letting  7,j  represent  the  intensity  output  value  of 
detector  pixel  ij  and  defining  the  x-direction  as  the  direction  moved  as  j is  incremented. 
The  x-directed  centroid  can  be  written  as 


Z 0'-^)  Z^/./ 

Z lA. 


[5.1] 


j=<I-[%\i=p-[%\ 

Likewise,  the  y-directed  centroid  is  found  by  reversing  the  order  of  the  summation,  and  is 
given  by 


‘>AV2\ 

Z 0-^)  ZA'.y 

oAVll  p\%\ 

Z llh, 

j=<l-[%\i=P-[%\ 


Figure  5-2  shows  the  algorithm  error  as  the  focal  spot  position  is  varied  across  a 
single  pixel  width  for  centroiding  with  odd  values  of  N.  These  curves  were  generated  by 
assuming  a sine-squared  diffraction  pattern  from  a square  lenslet  and  numerically 
integrating  the  intensity  values  to  find  the  detector  pixel  output  values.  The  position  of 
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the  spot  was  varied  by  plus  and  minus  one  half  of  a pixel  width  and  the  centroid, 
calculated  as  above,  was  subtracted  from  the  actual  position.  Note  that  the  error  is  zero 
when  the  spot  is  positioned  at  the  center  of  the  peak  reference  pixel  and  grows  to  a 
maximum  at  the  edge.  When  the  spot  position  crosses  over  into  the  adjacent  pixel,  the 
algorithm  error  will  swap  signs  resulting  in  a large  apparent  change  in  spot  position. 
Figure  5-3  shows  the  algorithm  error  as  the  focal  spot  position  is  varied  across  a single 
pixel  width  for  centroiding  with  even  values  of  N.  Note  that  the  error  is  zero  when  the 


Figure  5-2:  Odd  point  centroid  algorithm  errors 
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spot  is  positioned  at  the  boundary  between  two  pixels  and  is  a maximum  when  centered 
on  a pixel. 

In  addition  to  the  large  magnitude  of  the  algorithm  errors,  another  significant 
problem  with  the  centroiding  algorithms  is  the  discontinuity  in  the  error  as  a boundary  is 
approached.  The  odd-point  centroid  algorithm  uses  an  odd  number  of  detector  pixel  rows 
or  columns  consisting  of  the  row  or  column  with  the  peak  value  and  an  equal  number  on 
each  side  of  it.  As  seen  in  Figure  5-2,  the  odd-point  centroid  algorithm  has  a non-zero 
error  when  the  spot  is  centered  over  the  boundary  between  two  pixels.  The  calculated 
centroid  will  be  biased  toward  the  center  of  the  peak  pixel,  which  is  used  as  the  central 
pixel  in  the  calculation.  As  the  spot  shifts  slightly  beyond  this  position  such  that  the 


4,6,8  Even-Point  Centroid  Algorithm  Error 


Figure  5-3:  Even  point  centroid  algorithm  errors 
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adjacent  pixel  becomes  the  peak  reference  pixel,  the  calculated  position  becomes  biased 
toward  the  new  central  pixel.  This  results  in  a perceived  shift  equal  to  twice  the 
boundary  error  value.  The  even-point  centroid  algorithm  uses  an  even  number  of 
detector  rows  or  columns  containing  the  peak  row  or  column,  the  next  highest  row  or 
column,  and  an  equal  number  on  each  side  of  these  two.  It  therefore  has  zero  error  at  the 
boundary  between  two  pixels  and  non-zero  error  when  the  spot  is  centered  over  a pixel, 
as  seen  in  Figure  5-3.  When  the  spot  position  transitions  across  the  center  of  a pixel,  a 
row  or  column  on  the  side  from  which  it  is  moving  is  dropped  from  the  calculation  as  a 
new  row  or  column  on  the  side  toward  which  it  is  moving  is  added.  This  results  in  the 
discontinuity  shown  in  figure  5-3  at  zero  offset. 

At  first  glance,  the  types  and  magnitudes  of  algorithm  error  shown  and  discussed 
above  appear  to  be  prohibitively  large.  However,  since  these  errors  are  largely 
deterministic  in  nature,  it  is  possible  to  calibrate  out  much  of  the  error.  By  using 
centroiding  algorithms  which  are  single  valued  and  calculating  the  error  as  a function  of 
spot  position  as  done  in  Figures  5-2  and  5-3,  a look-up  table  can  be  generated  and  used  to 
correct  the  algorithm  error  calculations.  The  fidelity  of  the  algorithm  error  corrections 
depend  upon  how  closely  we  are  able  to  model  the  actual  diffraction  spot  profile  in  terms 
of  its  size  and  shape.  Figures  5-2  and  5-3  used  an  ideal  sine-squared  intensity  pattern, 
which  was  numerically  integrated  to  determine  the  detector  pixel  output  values.  Since 
the  actual  pattern  will  not  be  a perfect  sine-squared  pattern,  there  will  be  errors  in  the 
correction  data.  These  errors  will  be  most  apparent  at  the  transition  boundaries  where  the 
discontinuities  exist.  One  way  to  alleviate  this  problem  is  to  use  an  algorithm  which  does 
not  contain  these  discontinuities.  Figure  5-4  shows  the  algorithm  error  for  a 4 order 
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curve  fit  algorithm.  In  this  algorithm,  five  rows  or  columns  centered  on  the  peak  valued 
row  or  column  are  used  to  find  a fourth  order  polynomial  fit.  The  resulting  polynomial 
peak  value  and  its  position  can  be  found  by  determining  the  roots  of  its  derivative.  The 
most  important  feature  of  this  algorithm  is  the  fact  that  the  algorithm  error  is  zero  at  the 
center  position  as  well  as  at  the  boundaries.  Although  the  peak  magnitude  of  the  error  is 
larger  than  with  the  centroiding  algorithms,  the  lack  of  discontinuities  should  result  in  an 
algorithm  that  is  less  corrupted  by  the  errors  that  result  when  transitioning  between 
groups  of  detector  pixels. 


4th  Order  Curve  Fit  Algorithm  Error 


Spot  Position  (pixels) 


Figure  5-4:  4*’’  order  curve  fit  algorithm  errors 
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Noise  Induced  Error 

The  next  significant  source  of  spot  position  measurement  errors  is  the  noise  present 
on  the  detector  pixel  outputs.  Various  noise  sources  contribute  to  this,  including  shot 
noise,  read  noise,  dark  current  shot  noise,  and  non-uniformities  in  the  responsivity 
between  pixels.  Using  the  expressions  above,  for  the  various  position  measurement 
algorithms,  a detailed  noise  analysis  can  be  carried  out  to  derive  an  expression  for  the 
noise  in  the  position  measurement  due  to  the  noise  on  the  detector  pixel  outputs.  Figure 
5-5  shows  the  geometry  used  for  this  analysis.  The  detector  array  is  a 100%  fill  factor 
array  of  square  pixels  of  width  P.  The  pixels  are  ordered  with  row  index  i and  column 


Figure  5-5:  Detector  array  geometry 
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index  j.  An  N-point  centroid  calculation  will  be  performed  on  each  N by  N pixel  group, 
centered  on  the  pixel  containing  the  peak  value.  The  pixel  output  values  Sij,  represent  the 
responsivity-scaled  integrated  intensity  over  the  pixel  surface.  Each  centroid  will  be 
calculated  about  a reference  position  of  (0,0)  corresponding  to  the  center  of  the  pixel 
containing  the  peak  output  value.  For  the  following  analysis,  we  will  consider  the  x- 
directed  odd-point  centroid  calculation  corresponding  to  the  j direction.  We  can  more 
concisely  write  the  centroiding  equations  by  defining  the  column  sums  Sj,  as 

(V-l) 

Sj=  t S,,,.  [5.3] 

(V-1) 


We  can  now  write  the  centroid  as 

(Af-I) 

^ t iS, 

. (v-i) 

A = . [5.4] 

^ (V-l) 

t s, 

. (N-\) 


We  now  need  to  examine  the  various  noise  sources,  which  are  contained  in  the  Sj  terms. 
Shot  noise 

The  shot  noise,  which  is  due  to  the  discrete  nature  of  the  photons  comprising  the 
source,  is  dependent  upon  the  mean  number  of  photons  N ■ . If  we  lump  the  spectral 
responsivity  and  quantum  efficiency  terms  together  with  the  appropriate  scaling  factors, 
we  can  write  an  expression  for  detector  output  in  terms  of  the  number  of  incident  photons 
at  a fixed  wavelength  as  in  the  following  equation. 

^ Nj 


[5.5] 
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Since  the  RMS  shot  noise  level  in  the  source  is  equal  to  the  square  root  of  the  number  of 
photons  comprising  the  source 


The  read  noise  is  primarily  due  to  the  detector  array  output  amplifier  electronics. 
This  noise  source  will  have  a fixed  RMS  level  , which  is  camera  dependent. 

Dark  current  noise 

The  dark  current  noise  is  a shot  noise  term  due  to  the  dark  current  of  the  detector. 


The  dark  current  noise  term  is  typically  insignificant  when  compared  to  the  read  noise 
term  in  which  case  it  can  be  neglected. 

Responsivitv  non-uniformity 

The  non-uniformity  in  the  responsivity  01,  of  different  detector  pixels  can  be 
modeled  as  an  independent  noise  term.  If  we  define  the  uncertainty  of  the  responsivity  as 
a percentage  of  the  mean  responsivity  so  that, 


[5.6] 


Read  noise 


[5.7] 


[5.8] 


then  we  can  write  the  pixel  output  value  as 


Si  j=K0t'  ^^I{x,y)  dxdy  = K 01- 0.  J 


[5.9] 


Pixel  ij 


K is  a constant  and  I(x,y)  represents  the  intensity  over  the  pixel  surface  i,  j and  0ij  is  its 
integral  over  the  pixel  surface.  Since  ^is  the  only  random  variable  here  we  can  write 
the  following  relationships  for  mean  and  mean  square  values. 


We  can  now  write  the  variance  of  the  output  due  to  the  uncertainty  in  detector 
responsivity  as 


[5.11] 


This  can  be  arranged  to  a more  convenient  form  so  that  we  can  write 


[5.12] 


Typical  values  for  residual  responsivity  non-uniformity  are  0.2  to  2 percent,  and  as  a 
result  this  source  of  error  can  be  significant. 

Quantization  noise 

Since  the  detector  output  will  be  digitized,  it  will  contain  quantization  errors.  If 
the  discrete  output  step  size  is  A , and  we  represent  the  input  with  a uniform  probability 
density,  then  it  is  straightforward  to  show  that  the  standard  deviation  of  the  quantization 
noise  has  the  following  form. 


Given  the  source  and  nature  of  the  above  noise  sources  it  is  reasonable  to  model 
them  as  independent  noise  terms.  The  total  noise  variance  is  therefore  the  sum  of  the 
individual  variances. 


[5.13] 


[5.14] 


Note  that  for  general  purpose  non-cooled  detector  arrays  such  as  the  one  being 
used  in  this  work,  the  gain  factor  is  quite  small  ( ,j^0.0025  for  the  Kodak  Megaplus 
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1.4  camera).  Consequently,  the  shot  noise  term  is  insignificant  as  are  the  dark  current 
and  quantization  terms,  when  compared  to  the  read  noise  and  the  responsivity  non- 
uniformity. Also  note  that  while  the  read  noise  is  an  additive  term,  the  responsivity  non- 
uniformity noise  model  is  dependent  upon  the  signal  level. 

Centroiding  algorithm  noise  errors 

We  have  now  developed  an  expression  for  the  total  noise  variance  for  each  pixel 
output  value  or  column  sum  value.  We  will  now  examine  how  this  noise  term  propagates 
through  the  centroiding  equations  and  generates  noise  in  the  spot  position  measurement. 
Rewriting  equation  [5.4]  in  the  following  shorthand  notation  will  make  the  analysis  more 
tractable. 


Since  the  random  variable  Sj  is  contained  in  both  the  numerator  and  denominator  we 
cannot  write  a closed  form  solution  for  the  variance  of  the  spot  position  X,  without  further 
simplification.  We  will  expand  each  term  of  the  numerator  sum  in  a Taylor  series 

expansion  around  each  mean  value  and  drop  all  second  order  or  higher  terms. 


[5.15] 


[5.16] 


[5.17] 
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Now  if  we  find  the  variance  of  each  term  of  the  sum  over  n separately,  the  N 
variances  can  be  added  to  find  the  total  variance  in  X. 
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Now  the  total  variance  can  be  found  by  summing  the  N terms. 


f \ 

z^. 

V / y 


!■ 


,-|2 


n — 


nS, 


f \ 

V i 


[5.26] 


= 


( \ 

P 

f \ 

Y.S, 

VV ./  y 


[5.27] 


Note  that  equation  [5.27]  has  a number  of  variables  which  are  dependent  upon  the 
spot  position.  The  values  are  the  detector  pixel  column  sum  values,  which  obviously 
vary  as  the  spot  position  changes.  Likewise,  cr^  is  generally  dependent  upon  the  spot 


position  since  its  shot  noise  and  non-uniformity  terms  are  a function  of  as  described 
above.  The  dark  current  and  read  noise  components  of  are  constant  with  spot  position. 
Figure  5-6  shows  a plot  of  equation  [5.27]  normalized  to  P=l.  The  standard  deviation  of 
the  detector  outputs  were  set  by  a 20  dB  signal-to-noise  ratio  (SNR),  and  a null-to-null 
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Centroiding  Standard  Deviation  for  N=3,5,7 


Figure  5-6:  Odd-point  centroid  noise  standard  deviation  (SNR=20dB) 


spot  width  equal  to  2.9  times  the  pixel  width  was  used.  The  worst  case  standard  deviation 
occurs  at  the  edge  of  the  measurement  range  when  the  spot  is  straddling  two  adjacent 
pixels,  because  this  gives  the  largest  value  of  cr^  for  n^O.  Note  that  when  the  spot  is 

centered  on  the  central  pixel  the  cr^  value  is  maximum  for  that  pixel  but  it  gets  a 
weighting  value  of  n=0  and  thus  does  not  contribute  to  the  centroid  noise.  The  adjacent 
pixels  will  have  significantly  smaller  values  of  for  typical  spot  widths,  as  compared 
to  their  values  at  the  edge  of  the  measurement  range.  Therefore  the  minimum  value  for 
the  centroid  noise  standard  deviation  occurs  at  the  central  position.  Note  that  the  above 
example  represents  the  case  where  the  pixel  output  noise  is  a function  of  the  pixel  mean 
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output  value,  as  is  the  case  with  the  responsivity  non-uniformity.  If  a fixed  noise  level, 
such  as  read  noise,  is  dominant,  then  the  shape  of  the  curve  in  figure  5-6  will  be  more 
flat.  In  order  to  get  an  accurate  estimate  of  the  expected  noise  level  we  must  factor  in  the 
pertinent  camera  specifications  to  properly  model  the  detector  noise  sources.  This  will  be 
carried  out  later.  The  above  analysis  was  for  an  odd-point  centroid  algorithm.  By 
changing  the  limits  on  the  pixel  indices  appropriately,  the  same  equations  can  be  used  to 
find  the  standard  deviation  of  the  even-point  centroid  algorithm  noise. 

Curve  fitting  algorithm  noise  errors 

The  curve-fit  algorithms  discussed  previously  are  not  nearly  as  straightforward  as 
the  centroiding  algorithms  in  terms  of  the  noise  analysis.  The  basic  steps  in  these 
algorithms  consist  of  finding  a least  squares  polynomial  fit  to  the  data,  taking  its 
derivative,  and  finding  the  root  that  corresponds  to  the  central  peak.  In  order  to  examine 
how  the  detector  noise  affects  the  algorithm  accuracy,  it  is  much  more  straightforward  to 
model  the  detector  noise  and  perform  simulations  on  the  algorithm  in  order  to  generate 
empirical  data.  Figure  5-7  represents  the  results  for  this  type  of  simulation  noise  analysis. 
A fourth  order  curve  fit  algorithm  was  used.  The  noise  was  modeled  as  additive  gaussian 
noise  with  a standard  deviation  set  by  a fixed  SNR  according  to  the  detector  mean  output 
level.  For  each  spot  location,  noise  was  added  to  each  pixel  output  and  the  curve  fit 
algorithm  was  used  to  find  the  location  of  the  spot  peak.  1 000  manifestations  of  the 
noisy  outputs  were  taken  for  each  spot  position  and  the  standard  deviation  of  these  1000 
estimates  was  calculated.  Figure  5-7  shows  a plot  of  the  standard  deviation  of  the 
estimates  as  a function  of  spot  position  for  SNR  values  ranging  from  17  to  30  dB. 

Notice  that  these  curves  have  the  same  general  form  as  the  centroiding  algorithms  seen  in 
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Figure  5-7:  4*’’  order  curve  fit  noise  standard  deviation 


Figure  5-6,  with  a minimum  at  the  center  and  a maximum  at  the  edges. 

Misalignment  Errors 

Another  class  of  errors  in  the  Shack-Hartmann  sensor  is  the  misalignment 
between  the  lens  array  and  detector  array.  We  will  assume  that  if  there  is  any 
magnification  between  the  lens  array  and  detector  array,  then  it  is  perfect  and  can  be 
accounted  for  with  a simple  scaling  constant.  Imperfections  or  aberrations  due  to  the 
intermediate  magnification  optics  will  be  considered  in  the  next  section.  If  we  consider 
the  lens  array  to  be  fixed  within  the  system  coordinate  system,  then  misalignments  of  the 
detector  array  can  take  on  numerous  forms.  A detailed  analysis  of  all  possible 
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misalignments  was  carried  out  by  Pfund  et  al.  (1998).  In  this  work,  the  misalignments 
were  grouped  as  lateral  translation,  axial  translation,  rotation  around  a lateral  axis,  or 
rotation  about  the  axial  (optical)  axis.  In  our  work,  where  we  will  be  taking  a reference 
frame  of  data  prior  to  the  measurement,  the  lateral  translations  and  rotations  about  the 
optical  axis  do  not  come  into  play.  These  misalignments  result  in  fixed  offsets,  which  are 
present  and  unchanged  in  the  reference  and  test  frames.  Therefore,  when  the  reference 
frame  is  subtracted  from  the  test  frame  the  misalignment  errors  are  removed.  This  leaves 
axial  translation  and  rotation  around  a lateral  axis  for  further  consideration. 

Axial  translation  error 

Small  axial  misalignments  in  the  offset  between  the  lens  array  and  detector  array 
can  result  in  two  different  effects.  The  first  of  these  is  in  that  the  detector  array  does  not 
lie  in  the  focal  plane  of  the  lens  array  so  that  the  spot  diffraction  patterns  do  not  match  the 
models  used  to  develop  the  algorithm  error  correction  data.  However  because  the  lenslets 
used  in  these  applications  typically  have  very  large  f-number,  the  depth  of  field  is  quite 
large  and  reasonable  amounts  of  axial  misalignment  result  in  negligible  errors. 

Another  type  of  error,  which  can  result  from  axial  misalignment  is  an  incorrect 
scaling  between  spot  deflection  measurements  and  wavefront  tilt  angles.  As  discussed 
previously,  the  wavefront  tilt  angle  0,  is  related  to  the  spot  deflection  5,  by  the  lenslet 
focal  length  f. 

0 = - [5.28] 

/ 

This  assumes  however  that  the  spot  deflection  measurement  occurs  exactly  at  the 
focal  plane.  If  the  measurement  is  made  in  a plane  offset  from  the  focal  plane  by  some 


small  amount  of  axial  misalignment  e,  then  the  scaling  constant  is  changed  such  that. 
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e = ^ [5.29] 

f + £ 

In  our  work  we  will  eliminate  this  error  by  a calibration  procedure  which  measures  a 
known  set  of  linear  wavefront  tilts  using  the  shack-hartmann  sensor.  This  will  give  a 
direct  measure  of  the  scaling  constant  including  any  axial  misalignment  error. 

Rotation  about  a lateral  axis 

The  final  misalignment  type  to  be  considered  is  the  rotation  of  the  detector  array 
with  respect  to  the  lens  array  about  one  of  the  lateral  axes.  We  will  develop  a 
relationship  for  the  amount  of  error  caused  by  this  type  of  misalignment  in  order  to 
establish  a requirement  for  alignment  accuracy.  Referring  to  Figure  5-8,  which  shows  a 
misalignment  of  (5.  If  we  consider  the  chief  ray  for  a lenslet  at  some  offset  position  X, 
and  some  small  wavefront  tilt  angle  9,  then  the  spot  deflection  along  the  rotated  detector 


Lens  Array 


Rotated 
Detector  Array 

Unrotated 
Detector  Array 


Figure  5-8;  Lateral  rotation  misalignment  /3 
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array  Sr,  will  not  be  equal  to  the  expected  value  of f6  as  in  the  un-rotated  case.  Using  the 

geometry  of  Figure  5-8,  an  expression  for  4 can  de  developed. 

/-^tanyg  _(/-Ztan/?)sin^ 

sin  6*  cos{0-/3)  ''  cos{0-/^) 

The  error  in  the  spot  deflection  measurement  Sg , is  simply  the  difference  between  this 

and  the  ideal  value. 


^ {f-Xtanfi)sind^ 
^ ^ cos{0-f3) 


1- 


(l-^tan^) 
cos{6-  /3) 


s XOXd^nP 


[5.31] 


Once  the  system  requirements  for  spot  measurement  accuracy  and  detector  array  size  are 
known,  equation  [5.31]  can  be  used  to  establish  limits  on  lateral  rotation  misalignment. 
Optical  Aberration  Errors 

Any  optical  components  placed  between  the  lens  array  and  detector  array  will 
have  aberrations,  which  cause  deviations  in  the  ideal  spot  locations.  Most  of  these 
aberrations  will  not  have  a negative  effect  on  the  measurements  because  they  will  be 
present  and  unchanged  for  both  the  reference  and  test  frames.  Like  the  lateral  translation 
misalignments  errors,  they  will  be  subtracted  out  when  the  reference  frame  spot  locations 
are  subtracted  from  the  test  frame  to  find  the  spot  deviations.  The  primary  aberration  of 
concern  is  lateral  distortion,  which  results  in  a change  in  magnification  across  the  field  of 
the  measurement.  This  is  the  typical  definition  of  distortion  in  which  the  error  in 
refraction  power  is  only  a function  of  field  position.  This  type  of  distortion  will  cause 
different  lenslets  of  the  lens  array  to  have  different  spot  position  errors,  but  the  errors  for 
a given  lenslet  will  not  change  from  reference  to  test  frame  and  will  thus  be  eliminated. 

A more  obscure  form  of  distortion,  which  is  difficult  to  predict  or  to  quantify,  is  angle 
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dependent  distortion.  In  this  type  of  distortion  the  error  is  also  a function  of  the  field 
angle  of  arrival,  which  is  what  changes  between  reference  and  test  frame  and  is  the  focus 
of  our  measurement.  Careful  attention  to  this  type  of  distortion  will  be  in  order  if  any 
optical  components  are  to  be  used  between  the  lens  array  and  detector  array. 

Snot  Position  Measurement  Strategy 

Given  the  above  relationships,  we  can  now  turn  our  attention  to  the  spot  position 
measurement  problem.  This  problem  becomes  one  of  determining  the  size  and  number  of 
the  detector  pixels  to  use  in  a given  position  measurement  algorithm.  Although  the 
centroiding  algorithms  have  the  undesirable  discontinuity  characteristics  discussed  above, 
some  insight  is  gained  by  examining  their  error  characteristics  as  a function  of  spot  size 
and  number  of  pixels  used.  We  can  then  examine  the  curve  fit  peak  location  algorithms 
in  the  presence  of  detector  noise  and  residual  algorithm  error. 

Centroiding  Algorithms 

Figure  5-9  shows  a plot  of  the  worst  case  algorithm  and  noise  induced  error  as  a 
function  of  the  centroid  size  N,  and  the  null-to-null  spot  width  to  pixel  width  ratio.  The 
noise  level  of  each  pixel  is  set  by  a 20  dB  SNR  relative  to  the  pixel's  mean  level.  Since 
algorithm  error  correction  will  be  used  with  a goal  of  achieving  no  greater  than  5% 
residual  algorithm  error,  this  amount  of  algorithm  error  is  added  to  the  worst  case  noise 
induced  RMS  error  in  these  plots.  As  seen  in  the  figure,  the  error  decreases  with 
increasing  centroid  size,  N and  has  a minimum  for  a spot  size  to  pixel  width  ratio  of 
approximately  3.  The  noise  model  used  here  is  most  representative  of  the  responsivity 
non-uniformity  or  shot  noise  limited  case  since  the  noise  level  of  each  pixel  increases  as 
the  mean  level  of  that  pixel  increases.  The  read  noise  limited  case  would  have  a constant 
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Figure  5-9;  3-9  point  centroid  total  worst  case  RMS  error  vs.  spot  size 

(pixel  SNR=20  dB) 


noise  variance  for  every  pixel  regardless  of  the  pixel  values.  Figure  5-10  shows  the  same 
information  for  this  read  noise  limited  case  where  the  standard  deviation  of  the  noise  at 
every  pixel  is  set  at  20  dB  below  the  maximum  pixel  level.  Notice  that  this  case  indicates 
that  larger  spot  sizes  and  smaller  centroid  values  result  in  lower  errors,  although  there  is  a 
much  less  pronounced  dip  for  spot  width  to  pixel  width  ratios  near  3.5.  The  two  limiting 
cases  shown  in  figures  5-9  and  5-10  will  bound  the  actual  conditions  for  most  practical 
applications,  where  the  lower  mean  value  pixels  will  be  read  noise  limited  and  the  peak 
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(global  SNR=20  dB) 

pixel  values  will  have  increased  noise  due  to  responsivity  non-uniformity.  Without 
modeling  the  noise  characteristics  of  a specific  detector  array,  the  above  results  indicate 
that  a centroid  size  near  5 and  a spot  width  to  pixel  width  ratio  of  3 to  4 would  be  a good 
compromise. 

4^*^  Order  Curve-fit  Algorithm 

Since  the  curve  fit  algorithm  provides  a continuous  estimate  with  no 
discontinuities,  it  will  be  used  in  this  work.  The  4*'’  order  fit  is  selected  because  it  is 
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analagous  to  the  N=5  centroid  , which  had  good  error  characteristics,  and  because  an 
even  order  is  desired  to  model  the  diffraction  pattern  peak.  We  can  now  adjust  the 
detector  pixel  size  to  achieve  a good  balance  between  algorithm  error  and  noise  induced 
error.  Figure  5-1 1 shows  that  the  algorithm  error  is  decreased  as  the  spot  size  is 
increased  relative  to  the  detector  pixel  size.  The  peak  error  is  shown  here  since  it  is 


Worst  Case  Algorithm  Error  vs.  Spot  Size 


Figure  5-1 1 : 4*’’  order  curve  fit  worst  case  algorithm  error  vs.  spot  size 


possible  that  a measurement  could  be  made  with  a large  number  of  spot  positions  at  or 
near  the  peak  error  position.  If  algorithm  error  were  the  only  source  of  error,  then  we 
could  simply  make  the  spot  size  as  large  as  practical  to  minimize  the  noise.  However  the 
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errors  due  to  noise  must  be  considered  and  balanced  with  the  algorithm  errors.  Figure  5- 
12  shows  how  the  peak  noise  induced  error  is  affected  by  the  spot  size  when  the  noise 
level  is  at  a fixed  SNR  below  the  mean  detector  output  level  for  each  detector  pixel, 
which  corresponds  to  the  responsivity  non-uniformity  limited  case.  As  the  figure  shows, 


Worst  Case  RMS  Noise  Induced  Error  vs.  Spot  Size 


Figure  5-12:  order  curve  fit  peak  noise  induced  error  vs.  spot  size 

(pixel  SNR=17,  20,  23,  26,  30  dB) 


the  errors  increase  as  the  spot  size  is  increased  which  is  in  direct  contrast  with  the 
algorithm  error.  If  we  plot  the  total  error  by  adding  the  5 % residual  algorithm  error  to 
the  noise  induced  RMS  error,  then  the  curves  of  Figure  5-13  result.  The  curves  now  have 
a minimum  error  point,  which  changes  position  depending  upon  the  SNR.  The  minimum 
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Null-Null  Spot  Width/Pixel  Width 

Figure  5-13:  4*  order  curve  fit  total  worst  case  RMS  error  vs.  spot  size 
(pixel  SNR=17,  20,  23,  26,  30  dB) 


occurs  at  a null-to-null  spot  diameter  to  pixel  width  ratio  of  approximately  3.5  for  the  17 
dB  SNR  and  moves  out  to  around  5.0  for  high  SNR  values  approaching  30  dB.  Recall 
that  the  noise  induced  error  component  is  calculated  for  signal  dependent  noise  limited 
conditions.  If  the  detector  outputs  are  read  noise  limited,  then  the  noise  level  is  constant 
regardless  of  the  detector  output  level.  The  noise  induced  error  for  this  case  is  plotted  in 
Figure  5-14  in  which  the  SNR  values  set  the  standard  deviation  of  the  noise  for  all 
detector  pixels  equally.  The  SNR  sets  the  standard  deviation  of  the  read  noise  at  a fixed 
amount  below  the  peak  detector  pixel  output  value,  which  occurs  when  the  spot  is 
centered  on  a detector  pixel.  These  curves  also  have  a minimum  error  point,  which 
maintains  the  same  position  regardless  of  the  SNR.  The  minimum  occurs  at  a null-to-null 
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Figure  5-14:  4*  order  curve  fit  worst  case  RMS  noise  induced  error  vs.  spot  size 
(global  SNR=17,  20,  23,  26,  30  dB) 


spot  diameter  to  pixel  width  ratio  of  approximately  3.5.  If  we  plot  the  total  error  by 
adding  the  5 % residual  algorithm  error  to  the  read  noise  induced  error,  then  the  curves  of 
Figure  5-15  result.  The  curves  again  have  a minimum  error  point,  which  changes 
position  depending  upon  the  SNR.  The  minimum  occurs  at  a null-to-null  spot  diameter 
to  pixel  width  ratio  of  approximately  3.5  for  the  17  dB  SNR  and  moves  out  to  around  5.0 
for  high  SNR  values  approaching  30  dB.  The  minima  occur  at  nearly  the  same  points 
here  as  they  did  in  the  signal  dependent  noise  limited  case,  although  the  errors  are  slightly 
higher  in  this  read  noise  limited  case  for  equal  SNR  values.  The  curves  of  Figures  5-12 
through  5-15  demonstrate  the  effects  of  the  various  noise  sources  as  a function  of  the  spot 
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Total  Error  vs.  Spot  Size  (5%  Residual  Algorithm  Error) 


Figure  5-15:  4'^  order  curve  fit  total  worst  case  RMS  error  vs.  spot  size 
(global  SNR=17,  20,  23,  26,  30  dB) 


size.  However,  in  order  to  determine  the  optimum  spot  size  we  must  know  more  about 
the  actual  noise  characteristics  of  the  detector  array.  This  will  allow  us  to  model  all  of 
the  noise  sources,  and  determine  the  expected  errors. 

Expected  Performance  of  Kodak  Megaplus  1 .4  Camera 

In  this  work,  we  will  be  using  the  Kodak  Megaplus  1 .4  8-bit  camera  with  a 532 
nm  laser  as  the  source.  Table  5-1  outlines  the  pertinent  camera  characteristics. 

We  can  now  use  equation  [5.14]  to  set  the  noise  level  of  each  pixel  individually  and 
determine  the  total  error  as  a function  of  spot  size.  The  detector  array  noise 
characteristics  in  Table  5-1  were  used  to  model  the  noise  level  of  each  pixel  and  the  4*’’ 
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Table  5-1:  Kodak  Megaplus  1.4  camera  noise  characteristics 


Gain  Constant 

^1/392 

Read  Noise 

<0.5%  Full  Scale 

0--1.25 

Shot  Noise 

cr  = S 

Dark  Current  Noise 

Quantizer  Limited 

cr  =0.29 

Responsivity  non-Uniformity 

5 <0.35% 

a = SS 

order  curve  fit  algorithm  was  used.  The  worst  case  errors  were  then  calculated  as  a 
function  of  spot  size  and  plotted  in  Figure  5-16.  The  data  was  generated  using  the  same 
procedure  as  in  figures  5-12  and  5-14  except  that  the  noise  was  modeled  for  each  detector 
pixel  according  to  Table  5-1.  The  curve  in  figure  5-16  has  a minimum  error  at  a spot 
width  to  pixel  width  ratio  of  approximately  3.5  as  in  the  read  noise  limited  case  of  Figure 
5-14.  In  order  to  see  the  whole  picture,  we  must  also  consider  the  residual  algorithm 
error  in  conjunction  with  the  noise  error.  Assuming  that  the  algorithm  error  correction  is 
able  to  remove  all  but  5%  of  the  total  algorithm  error  shown  in  Figure  5-10,  the  noise 
induced  RMS  error  plus  worst  case  residual  algorithm  error  is  plotted  in  Figure  5-17. 

The  error  minimum  in  this  plot  has  moved  out  to  a spot  width  to  pixel  width  ratio  of 
approximately  4.0  where  the  total  error  is  approximately  1.3%  of  a pixel  width. 

In  the  previous  chapter,  it  was  shown  that  a 20  dB  dynamic  range  would  be 
required  in  the  measurement  of  wavefront  tilt  in  order  to  yield  sufficiently  accurate  phase 
reconstructions.  The  wavefront  tilt  measurement  corresponds  to  the  spot  position 
measurement  since  it  is  simply  a scaled  version  of  it.  The  above  analysis  indicates  that  a 
20  dB  dynamic  range  will  be  achieved  as  long  as  the  total  range  of  the  spot  position 
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Worst  Case  RMS  Noise  Induced  Enor  vs.  Spot  Size 


Figure  5-16:  4*’’  order  curve  fit  worst  case  RMS  noise  induced  error  vs.  spot  size 
(Kodak  Megaplus  1 .4  noise  model) 


measurement  is  greater  than  or  equal  to  1.3  pixels.  We  have  already  established  that  N=5 
because  the  4'^  order  curve  fit  algorithm  is  to  be  used.  If  a single  pixel  buffer  layer  is 
used  to  surround  the  N'  x N'  pixel  footprint  of  each  lenslet,  then  the  range  of  the  spot 
position  measurement  ± , will  be 


M-2 

2 


[5.32] 


Therefore  N'  > 5 results  in  a dynamic  range  which  exceeds  the  20  dB  requirement. 


Dynamic  Range  = 


AX 


-^  = 20.6  JR 
0.013 


[5.33] 
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Total  Error  vs.  Spot  Size  (5%  Residual  Algorithm  Error) 


Figure  5-17:  4*^  order  curve  fit  total  worst  case  RMS  error  vs.  spot  size 
(Kodak  Megaplus  1 .4  noise  model) 


Summary  and  Conclusions 

The  results  presented  in  this  chapter  describe  how  to  determine  the  spot  position 
measurement  accuracy  and  noise  characteristics.  The  results  presented  for  the  Kodak 
Megaplus  1 .4  camera  indicate  that  sufficient  dynamic  range  is  readily  achievable  using  a 
4*  order  curve  fit  algorithm.  This  algorithm  requires  a 5-by-5  group  of  pixels  centered 
on  the  peak  output  pixel  of  each  lenslet  footprint  group.  The  lenslet  footprint  should 
cover  an  M-by-M  group  of  detector  pixels  where  N'  is  at  least  5.  The  diffraction  pattern 
of  the  focussed  spot  should  have  a null-to-null  width  of  3.5  to  4.0  to  minimize  the  error. 
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These  results  will  be  used  in  the  next  chapter,  which  outlines  the  final  steps  in  the  Shack- 
Hartmann  system  design. 


CHAPTER  6 

SHACK-HARTMANN  SENSOR  DESIGN 
Problem  Statement 

In  the  previous  chapter,  we  examined  the  characteristics  of  the  detector  array  portion 
of  the  Shack-Hartmann  (S-H)  sensor.  We  will  now  use  this  knowledge  to  complete  the 
design  of  the  full  S-H  sensor,  which  includes  the  detector  array,  the  lens  array,  and  any 
required  intermediate  optical  components.  Given  a set  of  requirements  for  spatial 
resolution,  measurement  sensitivity,  dynamic  range,  and  test  region  size,  a suitable  lens 
array  and  detector  array  must  be  matched  to  meet  the  requirements.  The  resolution  of  the 
Hartmann  sensor  is  determined  by  the  sub-aperture  lenslet  spacing,  and  the  sensitivity  is 
determined  primarily  by  the  ability  to  measure  the  diffraction  spot  centroid  positions  on 
the  detector  array  in  the  presence  of  noise.  Unless  a custom  lens  array  is  to  be  used,  then 
the  lens  array  will  have  to  be  chosen  from  a limited  set  of  COTS  lens  arrays  so  that  it  is 
the  best  match  for  a given  detector  array.  Although  undesirable  in  terms  of  the  spot 
position  measurement  aberrations  it  produces,  a stage  of  imaging  optics  between  the  lens 
array  and  detector  array  provides  some  flexibility  in  matching  a lens  array  to  a detector 
array. 

Shack-Hartmann  Wavefront  Sensor  Design 
In  selecting  or  designing  the  lens  array,  the  key  parameters  are  lenslet  spacing 
(pitch),  lenslet  focal  length,  and  lens  array  size.  The  lens  array  samples  the  incoming 
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wave,  which  may  or  may  not  he  magnified  from  its  size  in  the  flow  test  region. 
Intermediate  magnification  or  demagnification  could  be  used  to  match  resolution 
requirements  to  available  lens  array  dimensions.  For  the  sake  of  completeness  we  will 
consider  the  general  case  shown  in  Figure  6-1,  where  there  is  magnification  AT  between 
the  test  region  and  the  lens  array,  as  well  as  magnification  M,  between  the  lens  array  and 


) 

) 

) 


Test  Region  Magnification  Lens  Magnification  Detector 

K Array  M Array 


Figure  6-1:  General  S-H  layout 

the  detector  array.  It  should  be  noted  that  the  magnification  between  the  test  region  and 
lens  array  must  preserve  the  phase  and  will  therefore  consist  of  a 4-f  double  lens 
arrangement. 

Key  Design  Relationships 

The  lens  array  pitch  A,  is  set  according  to  the  central  slice  theorem,  so  that  the 
resulting  reconstruction  has  the  required  spatial  resolution  (Kak  & Slaney  1988). 
Equation  [3.53]  can  be  modified  to  account  for  the  magnification  stage  K,  such  that 
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A < 


1.3823  K 


7Tf, 


[6.1] 


3dB 


In  this  application,  the  measurement  is  the  projected  phase  shift  through  the  test 
region  including  any  magnification  that  may  be  used  between  the  test  region  and  the  lens 
array.  Once  the  lens  array  pitch  has  been  established,  it  must  be  matched  to  the  detector 
array  and  the  total  lens  array  size  can  be  determined  to  yield  the  desired  test  region  size. 
Expected  wavefront  angles  (or  phase  gradients)  are  used  to  determine  the  required  focal 
length  for  the  required  sensitivity  and  dynamic  range.  Determination  of  the  required 
lenslet  focal  length  is  strongly  dependent  upon  the  spot  deflection  measurement  strategy. 
Following  the  approach  outlined  in  the  previous  chapter,  the  optimal  diffraction  spot  size 
^spot  determined  in  relation  to  the  detector  pixel  size,  P,  such  that 


=RP  . [6.2] 

where  R is  the  constant  of  proportionality  to  be  determined.  For  a square  lens  of  width  a 
and  focal  length  f the  diffraction  spot  null-to-null  diameter  for  light  at  wavelength  A can 
be  written  as 


^ spot 


lAf 


[6.3] 


Given  the  spot  size  and  detector  array  characteristics,  the  spot  position  measurement 
noise  variance  can  be  calculated  using  the  results  of  the  previous  chapter.  Keep  in  mind 
that  since  two  independent  measurements  are  being  subtracted  to  determine  the  spot 
deflection,  the  resulting  noise  variance  will  be  double  that  of  the  individual  spot  position 
measurements  discussed  in  Chapter  5.  The  spot  deflection  measurement  standard 
deviation  can  be  written  as  a fraction  of  the  detector  pixel  width  as 


cj^^aP  . 


[6.4] 
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This  will  set  the  minimum  measurement  sensitivity. 

The  maximum  measurement  range  is  determined  by  the  number  of  detector  pixels 
used  behind  each  lenslet.  In  order  to  avoid  overlap  between  adjacent  lenslets,  a boundary 
layer  of  at  least  one  detector  pixel  width  at  the  periphery  of  each  set  of  detector  pixels  is 
required.  Thus,  the  maximum  spot  deflection  measurement  5^^^  for  A^’ by  detector 
pixels  per  lenslet  can  be  written  as 
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Now  the  dynamic  range  requirement  D can  be  written  in  terms  of  the  number  of 
detector  pixels  per  lenslet  IST,  and  the  centroid  measurement  noise  factor  a , as 

d N'  -2 

D ^ ^ = i! — ± . [6.6] 

This  expression  can  then  be  used  to  determine  the  minimum  required  number  of  detector 
pixels  per  lenslet. 

N'>\2aD  + 2']  [6.7] 

Turning  our  attention  to  the  magnification  factor  M,  and  requiring  that  the  maximum 
wavefront  tilt  angle  yields  the  maximum  deflection  on  the  detector  array  results  in 
the  following  expression. 
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We  also  know  that  we  want  the  spot  diameter  on  the  detector  array  to  be  P/R,  so  that  we 
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Solving  equations  [6.8]  and  [6.9]  for / M and  equating  them  yields  an  expression  for  the 
lenslet  aperture  size. 


^max 

Now  by  matching  a total  of  Z lenslets  to  a total  of  Q detector  pixels,  an  expression  for  the 
total  magnification  can  be  derived. 


KM= 

XL{N'-2) 


[6.11] 


The  test  region  width  can  now  be  written  as 

^_La  ^LaM  ^ QP  ^XL{N'-2)  ^2] 

K KM  KM  R 0^,^ 

We  now  have  expressions  for  all  of  the  Shack-Hartmann  design  parameters.  Given  a set 
of  requirements  for  sensitivity,  dynamic  range,  and  spatial  resolution  the  wavefront 
sensor  can  now  be  designed  by  optimizing  the  various  magnifications,  centroid  sizes, 
lenslet,  and  detector  parameters  according  to  the  above  expressions. 

Design  Procedure 

The  basic  design  process  is  now  to  use  the  above  relationships  along  with  the 
design  requirements  and  component  constraints  to  find  the  optimum  lens  array  and 
magnification  stages  (K  and  M)  for  a given  detector  array  (camera).  The  optimum  design 
is  the  one  which  yields  the  maximum  dynamic  range,  while  meeting  or  exceeding  all 
design  requirements.  The  design  procedure  can  be  outlined  by  the  following  steps. 

Initial  iteration 


Using  the  known  requirements  and  the  characteristics  of  a given  spot  position 
measurement  algorithm,  determine  the  initial  estimates.  The  test  region  size  T,  is  a 
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Specified  requirement  and  the  sample  spacing  A,  can  be  determined  from  the  expected 
test  object  characteristics  using  equation  [6.1],  The  required  number  of  lenslets  L,  with 
fill  factor  , can  then  be  calculated  as 


The  results  of  the  previous  chapter  can  now  be  utilized  for  a given  spot  position 
measurement  algorithm.  The  optimum  spot  size  and  its  associated  RMS  error  can  be 
taken  from  plots  such  as  Figure  5-17  for  the  4^’  order  curve  fit  algorithm  or  Figures  5-9 
and  5-10  for  the  centroiding  algorithms.  The  optimum  RMS  error  factor  a,  along  with 
the  dynamic  range  requirement,  can  be  used  to  find  the  minimum  number  of  detector 
pixels  per  lenslet  N\  by  using  equation  [6.7].  Next,  the  total  magnification  can  be  found 
using  equation  [6.1 1].  At  this  point  we  have  initial  estimates  for  the  size  of  the  lens  array 
L,  the  number  of  detector  pixels  required  per  lenslet,  N',  and  the  total  magnification  KM. 
We  can  get  an  initial  estimate  for  K by  the  ratio  of  available  lens  array  sizes  to  test  region 
size.  For  example,  many  COTS  lens  arrays  are  available  in  25  mm  square  formats,  so  that 
a test  region  width  of  50  mm  establishes  a value  for  K of  0.5.  Given  this,  we  can  now 
find  the  lens  array  sub-aperture  size,  a.  We  can  also  estimate  the  lenslet  f-number  by 
requiring  that  the  maximum  wavefront  tilt  angle  results  in  a spot  deflection  near  the  edge 
of  the  lenslet  footprint.  For  example,  requiring  that  the  maximum  wavefront  tilt  causes  a 
spot  deflection  at  80%  of  the  footprint  half-width  results  in 


T dc^ 
A 


[6.13] 


= / 


0 


K 


max 


0.80  a 
2 


[6.14] 
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Lens  array  substitution  iterations 

Assuming  that  we  are  constrained  by  a particular  detector  array,  we  can  now 
begin  substituting  various  COTS  lens  arrays,  adjusting  K to  match  the  test  region  size  to 
the  lens  array  size,  and  adjusting  M to  match  the  lens  array  to  the  detector  array. 

Choosing  a lens  array  with  a size  L,  sub-aperture  diameter  a,  and  focal  length/,  which  are 
approximately  as  estimated  above,  we  can  determine  K and  M as  well  as  the  resulting 
spot  diameter,  dspot-  We  can  then  calculate  the  performance  characteristics  for  that  lens 
array  including  the  dynamic  range  D,  test  region  size  T,  and  maximum  wavefront  tilt 
9max-  This  can  be  done  with  varying  values  for  N'  and  the  best  combination  which  meets 
or  exceeds  the  requirement  will  be  chosen.  Note  that  equation  [6.6]  assumes  that  the  full 
range  of  spot  deflection  is  to  be  used.  Since  some  designs  may  not  make  use  of  the  full 
detector  sub-array,  we  need  to  write  the  dynamic  range  in  terms  of  the  maximum  utilized 
range.  The  dynamic  range  corresponding  to  this  maximum  wavefront  tilt  Omax,  will  be 

D = f ^ . [6.15] 

KaP 

A final  iteration  in  which  K is  adjusted  to  make  use  of  COTS  available  lenses  will  also  be 
required. 

Application  to  Design  Example 

Given  the  results  of  the  previous  chapters  along  with  the  design  procedure 
outlined  above,  we  can  now  design  the  S-H  wavefront  sensor  according  to  the 
requirements  of  the  design  example.  Since  this  design  will  be  used  to  perform  the 
experimental  phase  of  this  work,  it  will  be  constrained  by  the  practical  limits  imposed  by 
the  availability  of  COTS  lens  arrays  and  limited  CCD  camera  availability.  In  review,  the 
requirements  and  characteristics  of  the  design  example  are  as  follows. 
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Table  6-1:  Design  example  parameters 


Measurement 

Supersonic  Flow  over  an  Axi-symmetric 
Right  Circular  Cone 

Flow  Velocity 

Mach  1.5 

Shock  Angle 

60" 

Fluid 

Air  at  standard  conditions 

Test  Region  Size 

50mm  cube 

Reconstruction  resolution 

<lmm 

Optical  Wavelength 

532nm 

Maximum  Wavefront  Tilt  Angle 

0.2  ° (from  chapter  2) 

Projection  Data  Sample  Interval 

<0.44mm  (from  chapter  3) 

Spot  Position  Measurement  Dynamic 
Range 

>20  dB  (from  chapter  5) 

From  equation  [6.13],  and  assuming  a 100%  fill  factor,  the  required  number  of 
lenslets  will  be 


L = 


'T 

'50  X 1.0' 

A 

0.44 

= [113.6]  = 114 


[6.16] 


From  figure  5-17,  the  total  RMS  error  for  null-null  spot  widths  in  the  3 to  5 range 
will  be  approximately  0.015  times  the  pixel  width.  Multiplying  this  value  by  V2  and 
substituting  for  a in  equation  [6.7]  yields  a minimum  estimate  of  A^' . 

TV'  > [2  a D + 2]  = |’2xV2  X 0.015x100  + 2]  = [6.2]  = 7 [6.17] 

Given  a 1024  X 1024  detector  array  of  6.8//m  pixels  (Kodak  Megaplus  1.4  camera), 
equation  [6.1 1]  is  used,  along  with  an  estimate  for  R=A  (spot  diameter  equals  4 times  the 
pixel  width)  to  find  the  total  magnification  for  the  minimum  N'  of  7. 

QRPe^^^  1024x4x6.8x3.5x10'^ 


KM  = 


XRL{N' -2) 


0.32 


[6.18] 


0.532  X 114x5 
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Equation  [6.18]  gives  the  maximum  value  for  ATM  which  occurs  when  the  minimum  value 
of  N'  is  used.  In  order  to  match  the  lens  array  to  the  test  region  as  well  as  the  detector 
array,  a smaller  value  for  KM  may  be  required.  This  can  be  achieved  by  increasing  N' 
and/or  decreasing  R.  Exact  matching  of  the  50mm  test  region  to  a detector  array  of  1024, 
6.8mm  pixels  would  require  that. 


T K M = Q P=>  K M ^ 


Q± 

T 


1024x6.8 

50,000 


0.14  . 


[6.19] 


Equation  [6.19]  indicates  that  a larger  value  for  N'  and  smaller  value  of  will  be 
required.  An  estimate  of  the  lenslet  f-number  can  be  made  using  equation  [6.14], 
assuming  K=0.5. 


riu  0.40 0.40x0.5 
3.5x10-^ 

The  lenslet  aperture  diameter  for  K=0.5  would  be  as  follows, 

a = A K = 0.44  x 0.5  = 0.22mm 


[6.20] 


[6.21] 


and  the  focal  length  would  then  be, 

/ = «//#=  0.22x57  = 12.5mm  [6.22] 

Now  that  the  initial  iteration  is  complete,  we  can  begin  a search  for  lens  arrays, 
which  approximately  match  these  results.  A lens  array  with  approximately  1 14  X 114 
lenslets,  having  a lenslet  aperture  diameter  of  220//m,  and  focal  length  of  12.5mm  is  the 
target.  Once  we  find  those  lens  arrays  that  have  similar  specifications,  we  can  determine 
the  practical  values  for  K and  M to  match  the  lens  array  and  determine  the  expected 
performance  specifications  using  various  values  of  A^'.  For  our  design  example,  the 
following  lens  arrays  were  available  from  Adaptive  Optics  Associates. 
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Table  6.2  : Adaptive  Optics  Associates  COTS  lens  arrays 


ID 

Aperture  Diameter  (fwi) 

Focal  Length  (mm) 

Size 

1 

188 

8 

132  X 132 

2 

200 

6.5 

125  X 125 

3 

203 

7.8 

120  X 120 

4 

250 

18 

200  X 200 

As  an  example,  we  will  take  lens  array  1 through  the  steps  using  N'=l . In  each 
case,  we  will  use  all  available  lenslets  up  to  the  number  which  completely  fills  the 
detector  array  size  of  0=1024.  An  A' value  of  7 would  require  1024/7=146.3  lenslets  to 
fill  the  detector  array  completely.  Consequently,  the  maximum  available  number  of 
lenslets,  Z=132  will  be  used  with  lens  array  1.  Matching  the  lenslet  aperture  to  the 
detector  array  footprint  requires  that 


M = =0.253 

a 188 


[6.23] 


Now,  matching  the  50mm  test  region  to  the  lens  array  requires  that  the  first 
magnification  stage  have  magnification, 

K = — = = 0.496  , [6.24] 

T 50,000 


and  this  results  in  an  effective  sample  spacing  of 


a _ 0.188 
K ~ 0.496 


0.379  . 


[6.25] 
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The  spot  diameter  at  the  detector  array  can  now  be  found  using  equation  [6.3] 
along  with  the  calculated  value  for  M,  and  R can  be  determined  by  dividing  this  diameter 
by  the  pixel  width. 


d 


Spot 


2AfM  2x0.532x8000x0.253  ,,, 

= = 1 1.5  i/m 

a 188 


[6.26] 


R = 


11.5 

6.8 


= 1.7 


[6.27] 


The  results  of  chapter  5 can  now  be  used  to  determine  the  standard  deviation  of 
the  spot  position  measurement.  Using  the  data  in  Figure  5-17,  an  R value  of  1.7  results  in 


a standard  deviation  of  a = 0.023  V2  = 0.033  detector  pixels.  Equations  [6.6]  and  [6.15] 
can  now  be  used  to  calculate  the  resulting  dynamic  range. 


r.  / ^max  ^ ^OOO  X 3.5  X lO'"  X 0.253 

KaP  ~ 0.496  x 0.033  x 6.8  [6.28] 

Z)  = 63.6  = 18.1J5 


Finally,  we  can  calculate  the  maximum  wavefront  tilt  angle,  which  results  in  the 
maximum  useable  spot  deflection.  At  the  detector  array,  the  deflection  caused  by  a 
wavefront  tilt  of  Omax  after  both  stages  of  magnification  will  be 


5 = ^ . [6.29] 

max  ^ L -* 

We  can  now  determine  the  amount  of  wavefront  tilt  that  results  in  a spot  deflection  which 
is  one  detector  pixel  width  inside  the  lenslet  footprint. 


{N'-2)P 
K 2 


0. 


{N'-2)P  K 
2f  M 


[6.30] 


0 


max 


{N'-2)P  K 
2f  M 


5 X 6.8  X 0.496 
2 X 8000  X 0.253 


= A.llmrad  = 0.24° 


[6.31] 
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The  above  steps  were  applied  to  each  of  the  four  available  lens  arrays  using  N' 
values  of  7 to  12.  The  results  are  tabulated  in  Table  6-3.  The  only  case  in  which  the 
dynamic  range,  sample  spacing,  and  wavefront  tilt  specs  were  all  met  or  surpassed  is 
shown  in  bold  type.  In  selecting  the  best  choice  from  those  meeting  the  specifications,  it 
appears  that  lens  array  1 with  9 detector  pixels  per  lenslet  is  the  best  fit.  It  has  the  highest 
dynamic  range  while  meeting  the  sample  spacing  specification,  and  exceeding  the 
maximum  wavefront  tilt  slightly.  Also  important  is  the  fact  that  both  stages  of 
magnification  are  reasonable  in  magnitude,  and  somewhat  evenly  distributed. 

We  now  need  to  perform  a final  iteration,  while  taking  the  practical  limitations  of 
the  magnification  optics  implementation  into  account.  The  final  stage  of  magnification  is 
a simple  image  scaling  problem,  which  could  be  easily  implemented  using  a variable 
focus  camera  lens,  assuming  that  the  lens  has  enough  range  of  focus.  In  our  case  the 
Nikor  50mm  f/1.4  lens  does  not  so  that  an  additional  lens  will  be  required  between  the 
lens  array  and  the  camera  lens  in  order  to  image  the  array  of  focussed  spots  at  the  desired 
value  of  M=0.326.  A 160mm  focal  length  doublet  placed  approximately  121  mm  from 
the  lens  array  focal  plane  and  1 83mm  from  the  camera  lens  provides  the  desired 
magnification  with  this  camera  lens.  The  initial  magnification  stage,  K on  the  other  hand, 
must  be  a 4-f  system  in  order  to  preserve  the  wavefront  phase.  In  order  to  avoid  the 
requirement  for  custom  lenses,  it  is  desireable  to  make  use  of  COTS  available  lenses. 
Cemented  doublets  will  be  used  here  to  minimize  spherical  aberrations.  The 
magnification  of  a double  lens,  4-f  system  is  simply  the  ratio  of  the  focal  lengths,  and  the 
desired  magnification  factor  is  0.425  as  seen  in  Table  6.3.  Since  a limited  number  of 
focal  lengths  are  commercially  available  as  COTS  lenses,  a "best  fit"  design  will  be 
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required.  In  this  case,  the  best  match  is  achieved  with  a 500mm  and  200mm  lens  pair, 
which  gives  a magnification  factor  of  0.4.  Since  the  magnification  power  has  been 
decreased,  a larger  test  region  can  now  be  accommodated  using  the  same  number  of 
lenslets.  However  some  of  the  performance  characteristics  will  also  change,  and  a final 


Table  6-3:  COTS  lens  array  substitution  results 


ID 

N' 

L 

M 

R 

a 

K 

A 

D 

^max 

1 

7 

132 

0.253 

1.7 

0.033 

0.496 

0.38 

18.1 

0.24 

1 

8 

128 

0.289 

1.9 

0.030 

0.481 

0.39 

19.2 

0.24 

1 

9 

113 

0.326 

2.2 

0.028 

0.425 

0.44 

20.5 

0.22 

1 

10 

102 

0.362 

2.4 

0.025 

0.384 

0.49 

21.8 

0.21 

1 

11 

93 

0.398 

2.6 

0.024 

0.350 

0.54 

22.7 

0.19 

1 

12 

85 

0.434 

2.9 

0.021 

0.320 

0.59 

23.7 

0.18 

2 

7 

125 

0.238 

1.2 

0.037 

0.500 

0.40 

16.4 

0.31 

2 

8 

125 

0.272 

1.4 

0.035 

0.500 

0.40 

17.1 

0.33 

2 

9 

113 

0.306 

1.6 

0.034 

0.452 

0.44 

18.2 

0.31 

2 

10 

102 

0.340 

1.7 

0.033 

0.408 

0.49 

19.3 

0.29 

2 

11 

93 

0.374 

1.9 

0.030 

0.372 

0.54 

20.5 

0.27 

2 

12 

85 

0.408 

2.1 

0.028 

0.340 

0.59 

21.5 

0.25 

3 

7 

120 

0.234 

1.4 

0.035 

0.487 

0.42 

17.4 

0.26 

3 

8 

120 

0.268 

1.6 

0.034 

0.487 

0.42 

18.1 

0.27 

3 

9 

113 

0.301 

1.8 

0.031 

0.459 

0.44 

19.3 

0.27 

3 

10 

102 

0.335 

2.0 

0.030 

0.414 

0.49 

20.4 

0.25 

3 

11 

93 

0.368 

2.2 

0.028 

0.378 

0.54 

21.4 

0.23 

3 

12 

85 

0.402 

2.4 

0.025 

0.345 

0.59 

22.6 

0.21 

4 

7 

146 

0.190 

2.1 

0.028 

0.730 

0.34 

19.3 

0.21 

4 

8 

128 

0.218 

2.5 

0.025 

0.640 

0.39 

20.7 

0.19 

4 

9 

113 

0.245 

2.8 

0.023 

0.565 

0.44 

21.9 

0.17 

4 

10 

102 

0.272 

3.1 

0.021 

0.510 

0.49 

22.8 

0.16 

4 

11 

93 

0.299 

3.4 

0.020 

0.465 

0.54 

23.6 

0.15 

4 

12 

85 

0.326 

3.7 

0.018 

0.425 

0.59 

24.3 

0.14 

iteration  is  required  to  ensure  that  the  specifications  are  still  met.  The  new  test  region 
size  can  be  found  simply  as 
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L a 

Ic 


113x0.188 

0.4 


SIXmm  . 


[6.32] 


This  results  in  an  effective  sample  spacing  of , 


K 0.4 


[6.33] 


The  dynamic  range  can  be  calculated  as  before,  substituting  the  new  K value. 


^ , 8000  X 3.5  X 10-  X 0.326  ^ ^ ^ 


KaP 


0.4  X 0.028  X 6.8 


[6.34] 


Finally,  the  maximum  measureable  wavefront  tilt  angle  is 


= 


{N'-2)P  K _ 7 X 6.8  X 0.4 

2/  M 


2 X 8000  X 0.326 


= 3.65mrad  = 0.21° 


[6.35] 


Figure  6-2  shows  the  Shack-Hartmann  sensor  implementation  as  designed.  Note 
that  the  effective  sampling  interval  of  0.47mm  is  slightly  higher  than  the  0.44  mm 
requirement  calculated  earlier,  because  of  the  need  to  use  COTS  lenses.  Equation  [3.53] 


f= 500mm 

t 

53.1mm 

i 


Kodak  Megaplus  1.4 
21.2mm  with  f=50mm  lens 

1024X1 024  pixels 

f=200mm  1 f=  160mm  


500mm 


500mm 


X - > < 8mm 

200mm  200mm 

Lens  Array 
f=8mm 
L=113 
a=188um 


Figure  6-2:  Shack  Flartmann  sensor  implementation 
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can  be  used  again  to  show  that  the  spatial  resolution  will  now  be  1.07  mm,  which  is  out  of 
spec  by  7%.  If  the  spatial  resolution  requirement  is  critical,  then  a different  design 
should  be  considered,  possibly  including  a custom  lens  to  achieve  the  desired 
magnification,  K=0.425.  For  our  work,  the  1 mm  spatial  resolution  spec  was  somewhat 
arbitrary  and  so  the  design  of  figure  6-2  will  be  used. 

Summary  and  Conclusions 

This  chapter  has  outlined  the  Shack-Hartmann  design  process  as  applied  to  the 
aerodynamic  flow  density  measurement  problem.  The  design  procedure  was  tailored  to 
make  use  of  a particular  camera  as  well  as  non-custom  optics  and  lens  arrays.  In 
applications  where  these  constraints  are  not  present  or  where  a particular  design 
specification  requires  custom  components,  the  same  basic  design  relationships  can  be 
used.  The  above  procedure  was  developed  for  the  general  case  where  optical 
magnification  stages  were  used  to  match  the  test  region  to  the  lens  array  and  the  lens 
array  to  the  detector  array.  Because  of  the  potential  for  aberration  errors,  these 
components  would  ideally  be  eliminated  by  the  use  of  matching  lens  array,  detector  array 
pairs.  The  same  design  procedure  can  be  used  for  cases  in  which  they  are  not  present,  by 
setting  the  magnification  factors  to  one. 
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CHAPTER  7 

MEASUREMENTS  / RESULTS 
Problem  Statement 

Given  the  prototype  Shack-Hartmann  sensor  designed  from  the  previous  chapter, 
the  next  task  is  to  devise  a method  for  testing  its  measurement  characteristics.  The 
characteristics  of  critical  interest  include  measurement  accuracy,  sensitivity,  dynamic 
range,  and  repeatability.  The  simplest  and  most  effective  way  to  test  for  these 
performance  measures  is  to  utilize  a linear  wavefront  phase  or  constant  wavefront  tilt 
angle.  In  this  type  of  test,  the  measurements  from  each  lenslet  of  the  array  should  be  the 
same.  Assuming  that  a linear  wavefront  phase  is  presented  to  the  sensor,  any  deviations 
from  lenslet  to  lenslet  will  be  due  to  the  sensor  itself  It  will  then  be  possible  to 
determine  the  accuracy,  sensitivity  and  dynamic  range  by  varying  the  wavefront  tilt  angle 
over  the  full  dynamic  range  of  the  sensor.  An  independent  measure  of  the  wavefront  tilt 
will  be  required  if  the  accuracy  is  to  be  determined,  and  this  could  be  accomplished  in 
several  ways  as  discussed  later.  The  repeatability  can  be  determined  by  setting  the 
wavefront  tilt  angle,  taking  multiple  measurements  over  an  extended  time  period,  and 
comparing  them.  If  there  are  any  type  of  short  or  long  term  stability  issues,  they  will 
become  apparent  with  this  test. 

As  a demonstration  of  the  sensor's  capabilities  in  a dynamic  measurement 
application,  it  was  desired  to  perform  a measurement  such  as  the  supersonic  conical  flow 
for  which  the  sensor  was  designed.  However,  due  to  the  equipment  limitations  under 
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which  this  work  was  performed,  this  was  not  possible.  The  frame  grabber  hardware  and 
software  interfaces  available  for  the  Kodak  Megaplus  1 .4i  camera  did  not  have  external 
trigger  capability  and  a major  upgrade  was  required  to  replace  the  obsolete  frame  grabber 
board.  Because  of  time  and  budget  constraints  an  alternative  dynamic  demonstration 
measurement  was  chosen.  As  a simple  demonstration  of  a flow  pattern,  which  is  familiar 
to  the  aerodynamics  community,  a free  air  jet  nozzle  flow  measurement  will  be  made. 

This  measurement  can  be  made  by  setting  up  a quasi-steady  state  flow  condition,  thus 
alleviating  the  need  for  the  external  trigger  capability.  Although  the  wavefront  sensor 
was  not  designed  for  this  measurement  application,  it  is  hoped  that  its  inclusion  will 
provide  some  insight  into  the  operation  of  the  sensor  with  an  already  familiar  flow 
pattern. 

Static  Testing 

Experimental  Setup 

There  are  two  options  for  setting  the  wavefront  tilt  angle  as  discussed  above.  The 
option  chosen  for  this  work  is  to  use  a simple  adjustable  tilt,  folding  mirror  ahead  of  the 
wavefront  sensor.  A simple  micrometer  adjustment  allows  a continuous  adjustment  of 
the  wavefront  tilt  in  one  plane  of  the  wavefront  sensor.  An  independent  measurement  of 
tilt  angle  is  still  required  and  is  accomplished  by  beam-splitting  a portion  of  the  tilted 
beam  and  using  a single  lens  along  with  a linear  detector  array.  The  single  lens  and  linear 
detector  array  is  essentially  a single  element,  high  resolution,  one-dimensional  Shack- 
Hartmann  sensor.  By  focussing  the  beam  onto  the  linear  detector  array,  the  spot 
deflection  can  be  measured  to  determine  the  wavefront  tilt  angle. 
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The  experimental  set-up  used  in  this  work  is  shown  in  Figure  7-1 . The  4/ 
magnification  stage  between  the  test  region  and  the  lens  array  is  used  to  preserve  the 
phase  of  the  test  beam.  It  is  set  up  so  that  the  image  plane  is  located  at  the  lens  array  and 
the  object  plane  is  located  at  the  tilt  axis  of  the  adjustable  mirror.  This  is  critical  since  it 
insures  that  there  is  minimal  translation  of  the  beam  at  the  lens  array  as  the  tilt  angle  is 
adjusted.  A linear  detector  array  and  a 1500w/n  focal  length  lens  were  used  as  an 
independent  means  of  measuring  the  wavefront  tilt.  The  linear  detector  array  had  a Xh/dm 


Figure  7-1:  Experimental  setup  for  static  tests 
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pixel  pitch,  such  that  a single  pixel  of  measurement  accuracy  corresponded  to  0.009  mrad 
of  wavefront  tilt  measurement  accuracy.  This  is  a factor  of  4 lower  than  the  wavefront 
sensor  sensitivity.  All  lenses  shown  other  than  the  camera  lens  were  COTS  achromats. 
The  camera  used  a Nikor  50mm  f\A  lens. 

Another  option  for  setting  the  wavefront  tilt  angles  would  be  to  mount  the  lens 
array  and  all  optics  behind  it,  including  the  camera,  in  a rigid  housing  allowing  the 
wavefront  tilt  angle  to  be  set  by  tilting  this  sensor  assembly,  instead  of  the  wavefront. 

The  axis  of  rotation  would  be  in  the  plane  of  the  lens  array.  A housing  of  this  sort  would 
be  required  in  an  operational  system  anyway,  and  the  optics  used  to  form  the  test  beam 
become  much  simpler  to  design  and  build.  However,  due  to  the  tight  mechanical 
tolerances  required  to  fabricate  this  part,  and  the  limitations  of  available  machining 
capabilities  this  path  was  not  chosen. 

Test  Procedure 

The  wavefront  tilt  in  a single  plane  was  adjusted  across  the  entire  dynamic  range 
of  the  wavefront  sensor  using  the  tilt  mirror  and  the  linear  detector  array  and  camera  data 
was  recorded.  A reference  frame  of  data  was  recorded  by  finding  the  focussed  spot 
centroid  positions  with  the  wavefront  tilt  angle  set  to  approximately  zero.  These 
reference  centroid  positions  were  used  as  the  origins  for  all  subsequent  measurements. 
Processing  of  each  of  these  measurements  resulted  in  113^  wavefront  tilt  measurements 
in  both  directions  evenly  spaced  across  the  53.5  mm  square  test  region.  By  subtracting 
the  reference  centroid  positions  from  the  shifted  centroid  positions,  the  shifts  in  spot 
positions  were  calculated.  The  mean  value  and  standard  deviation  of  the  wavefront  tilts 
were  calculated  and  compared  to  the  wavefront  tilt  measured  by  the  linear  detector  array. 
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Measurement  Results 

Figure  7-2  shows  the  measured  wavefront  tilt  as  a function  of  the  actual  tilt,  as 
measured  by  the  linear  detector  array.  Figure  7-3  shows  the  average  error  in  measured 
tilt,  while  Figure  7-4  shows  that  the  standard  deviation  is  out  of  spec  for  most  wavefront 
tilt  angles.  Figure  7-5  is  a plot  of  all  wavefront  tilt  measurements  and  shows  that  there  is 
a distortion  error,  which  increases  significantly  toward  the  edges  of  the  wavefront  sensor. 
Further  analysis  by  ZEMAX  optical  design  software  has  indicated  that  this  is  a field 
angle  dependent  distortion  caused  by  the  optics  between  the  lens  array  and  the  wavefront 
sensor.  In  order  to  verify  this  experimentally,  the  wavefront  tilt  angle  was  set  to  zero  and 
the  lens  array  was  translated  normally  to  the  optical  axis  across  the  full  range  of  the 
measurement,  using  the  setup  of  Figure  7-6.  This  allowed  a measurement  of  spot 
position  to  be  made  without  changing  the  wavefront  tilt  angle.  Figure  7-7  shows  that  the 
standard  deviations  of  these  measurements  were  well  within  spec  for  the  full  dynamic 
range.  Figure  7-8  shows  that  the  angle  dependent  distortion  is  no  longer  present.  The 
only  difference  between  this  measurement  and  the  actual  measurement  is  that  the  angle  of 
the  bundle  of  rays  associated  with  each  lenslet  is  always  parallel  to  the  optical  axis  as 
opposed  to  tracking  the  wavefront  tilt  angle.  This  is  represented  in  Figure  7-9.  Although 
the  rays  associated  with  a lenslet  at  the  edge  of  the  lens  array  still  fall  within  the  entrance 
pupil  of  the  relay  optics,  the  distortion  which  is  more  a function  of  ray  angle  than  field 
position,  is  significant  as  the  edges  of  the  wavefront  sensor  are  approached.  Further 
improvement  of  the  optics  was  not  deemed  to  be  a viable  option,  since  the  required 
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Figure  7-2  ; Ideal  and  measured  wavefront  tilt  vs  actual  tilt 


Wavefront  Tilt  (mrad) 


Figure  7-3  : Measured  wavefront  tilt  error 
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Figure  7-4:  Standard  deviation  of  wavefront  tilt  measurements 
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Figure  7-5:  Spot  deviation  measurements  for  -3.1  mrad  wavefront  tilt 
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Figure  7-6  : Lens  array  translation  test  setup 


distortion  measure  at  the  edge  of  the  wavefront  sensor  would  be  on  the  order  of  0.007%. 
As  a result,  one  of  the  key  conclusions  of  this  work  is  the  recommendation  to  eliminate 
any  optics  between  the  lens  array  and  the  detector  array,  which  will  normally  mean  using 
a custom  lens  array.  In  order  to  make  some  use  of  the  current  set-up,  a subset  consisting 
of  the  35  X 35  central  measurements  were  used.  Figure  7-10  shows  the  measured  vs. 
actual  wavefront  tilt  angles  and  the  average  error  is  plotted  in  Figure  7-11.  Figure  7-12 
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Figure  7-7  ; Lens  array  translation  test  results 


Figure  7-8:  Lens  array  translation  test  results-spot  deviation  measurements 
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Figure  7-9  : Spot  translation  measurements  due  to  wavefront  tilt  and 
translation  of  lens  array 


shows  that  the  standard  deviations  for  these  measurements  are  within  the  specified 
requirement.  Figure  7-12  shows  that  the  angle  dependent  distortion  is  still  a significant 
source  of  error.  In  order  to  determine  how  stable  the  measurements  are  in  terms  of 
repeatability,  a set  of  time  delayed  measurements  were  performed.  The  wavefront  tilt 
angle  was  set  to  approximately  zero  and  a reference  frame  of  data  was  taken.  Additional 
frames  of  data  were  taken  at  time  delays  of  15  seconds,  2 minutes,  10  minutes,  and  60 
minutes.  The  spot  centroid  deviations  and  corresponding  wavefront  tilt  angles  from  the 
reference  frame  were  calculated  over  the  full  1 13  by  113  array.  Any  time  instabilities  in 
the  measurement  process  should  be  evident  in  the  means  and  standard  deviations  of  the 
delayed  measurements.  Figure  7-14  shows  the  mean  value  of  the  measured  wavefront  tilt 
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Figure  7-10:  Measured  vs  actual  tilt  for  35  X 35  data  subset 
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Figure  7-11:  Measured  tilt  error  for  35  X 35  data  subset 


Measured  Tilt  Standard  Deviation  (mrad) 


142 


-4-3-2-10  1 2 3 4 

Wavefront  Tilt  (mrad) 

Figure  7-12:  Standard  deviation  of  measured  tilt  for  35  X 35  data  subset 
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Figure  7-13:  Spot  deviation  measurements  for  -3.1  mrad  wavefront  tilt 

(35  X35  data  subset) 
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Figure  7-14  : Time  delayed  wavefront  tilt  means  (x  and  y directions) 
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Figure  7-15:  Time  delayed  wavefront  tilt  standard  deviations  (x  and  y directions) 
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and  Figure  7-15  shows  that  the  standard  deviation  for  the  delayed  measurements  are  all 
well  below  the  minimum  sensitivity  spec  of  0.035  mrad.  As  seen  in  Figure  7-14,  the 
mean  tilt  values  begin  to  exceed  the  0.035  mrad  sensitivity  spec  for  time  delays  over 
about  10  minutes.  This  should  not  be  surprising  since  no  effort  was  made  to  maintain  a 
stable  temperature  environment.  This  amount  of  wavefront  tilt  shift  corresponds  to 
approximately  0.3  microns  of  shift  in  spot  position  on  the  detector  array.  Considering  the 
nature  of  the  benchtop  setup  used,  changes  in  room  temperature  over  such  an  interval 
could  easily  cause  this  type  of  change.  Note  that  in  practice,  a reference  frame  of  data 
would  typically  be  taken  seconds  before  the  test  frame.  This  fact,  coupled  with  the  added 
mechanical  benefits  of  a rigid  body  mount  should  alleviate  any  concerns  over  the  stability 
and  repeatability  of  the  measurement. 

Dynamic  Testing 

Experimental  Setup 

In  order  to  demonstrate  the  use  of  the  system  in  a somewhat  familiar 
measurement  application,  a free  air  jet  from  a converging-diverging  nozzle  was  used.  A 
nozzle  was  designed  to  deliver  a mach  1.8  exit  velocity.  The  nozzle  design  is  shown  in 
Figure  7-16.  The  initial  plan  was  to  operate  the  nozzle  at  the  design  pressure  to  establish 
a steady  state  flow  at  the  design  conditions.  However,  a lack  of  compressed  air  capacity 
did  not  allow  for  this  to  be  accomplished.  Instead,  a tank  was  filled  to  80  psig  and 
allowed  to  empty  through  the  nozzle  upon  opening  an  exit  valve.  This  setup  is  shown  in 
Figure  7-17.  With  the  nozzle  in  the  test  region,  test  data  was  taken  as  the  pressure 
dropped  slowly  from  80  psig  to  10  psig. 
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Figure  7-16:  Convergiong-diverging  nozzle  design 


Figure  7-17  : Free  air  jet  nozzle  and  tank 


Test  Procedure 

The  first  step  in  the  measurement  procedure  was  to  take  a reference  frame  of  data 
with  the  nozzle  absent  from  the  test  region.  The  tank  and  nozzle  were  then  moved  into 
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position  so  that  the  nozzle  exit  was  just  visible  in  the  test  region  as  shown  in  Figure  7-18. 
After  filling  the  tank  to  80  psig,  the  exit  valve  was  fully  opened  and  the  camera  was 
manually  triggered  as  the  pressure  dropped  through  the  desired  pressure  reading.  Due  to 
the  rapid  initial  pressure  drop,  the  highest  tank  pressure  tested  was  60  psig.  Additional 


Figure  7-18  : Free  air  jet  test  setup 


measurements  at  10  psig  increments  from  60  to  10  psig  were  performed.  The  camera 
frame  time  was  set  to  its  minimum  value  of  10  msec. The  camera  data  was  processed  in 
the  same  manner  as  used  in  the  static  tests.  The  spot  deflections  were  calculated  and  used 
to  reconstruct  the  phase  map.  The  piecewise  zonal  reconstruction  discussed  in  Chapter  4 
was  used. 

After  all  cases  were  tested  using  the  wavefront  sensor,  dark  central  ground  and 
Schlieren  images  were  made  at  each  pressure  setting  for  comparison  purposes.  The  dark 
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central  ground  images  were  made  by  simply  removing  the  lens  array  and  placing  a 50 
micron  reflective  spot  at  the  focus  of  the  4f  demagnification  stage.  The  schlieren  images 
were  similarly  made  by  placing  a knife-edge  at  this  same  location. 

Measurement  Results 

In  the  following  figures,  the  reconstructed  phase  is  plotted  for  each  of  the  tested 
tank  pressure  readings.  Each  plot  of  reconstructed  phase  is  followed  by  the  dark  central 
ground  and  schlieren  images  for  that  pressure  setting.  It  should  be  noted  that  these  three 
measurements  were  not  simultaneously  performed  and  that  the  camera  was  manually 
triggered  as  the  tank  pressure  dropped  through  the  desired  value.  Consequently,  it  is 
likely  that  the  actual  tank  pressure  could  have  varied  by  as  much  as  2 psi  between 
measurements. 

As  seen  in  the  figures,  the  flow  for  tank  pressures  of  30  to  60  psig  is 
overexpanded  so  that  the  exit  Mach  number  will  be  at  approximately  1.8.  For  tank 
pressures  below  30  psig,  no  shocks  are  visible  indicating  subsonic  flow  conditions.  Mach 
numbers  for  these  conditions  will  be  bounded  by  an  upper  limit  of  0.62.  These 
approximations  are  based  on  one  dimensional  isentropic  and  normal  shock  flow 
assumptions. 

The  measurements  performed  above  represent  a single  projection  of  the  full  set  of 
data,  which  would  be  required  to  perform  the  tomographic  reconstruction.  However,  by 
assuming  that  the  flow  is  axially  symmetric,  it  is  possible  to  perform  this  final  step.  It 
should  be  noted  though,  that  the  phase  shift  data  above  is  only  somewhat  axially 
symmetric  and  that  the  noise  is  no  longer  uncorrelated  between  projections,  so  that  the 
tomographic  reconstructions  will  contain  significant  errors.  The  resulting  density  ratios 
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Figure  7-19  : Reconstructed  phase  shift  (a)  and  iso-phase  contours  (b)  tank  pressure  of  60  psig 
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Figure  7-20  : Dark  central  ground  (a)  and  Schlieren  (b)  images  for  tank  pressure  of  60  psig 
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Figure  7-21  : Reconstructed  phase  shift  (a)  and  iso-phase  contours  (b)  tank  pressure  of  50  psig 
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Figure  7-22  : Dark  central  ground  (a)  and  Schlieren  (b)  images  for  tank  pressure  of  50  psig 
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Figure  7-23  : Reconstructed  phase  shift  (a)  and  iso-phase  contours  (b)  tank  pressure  of  40  psig 
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Figure  7-24  : Dark  central  ground  (a)  and  Schlieren  (b)  images  for  tank  pressure  of  40  psig 
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Figure  7-25  : Reconstructed  phase  shift  (a)  and  iso-phase  contours  (b)  tank  pressure  of  30  psig 
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Figure  7-26  : Dark  central  ground  (a)  and  Schlieren  (b)  images  for  tank  pressure  of  30  psig 
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Figure  7-27  : Reconstructed  phase  shift  (a)  and  iso-phase  contours  (b)  tank  pressure  of  20  psig 
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Figure  7-28  : Dark  central  ground  (a)  and  Schlieren  (b)  images  for  tank  pressure  of  20  psig 
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Figure  7-29  : Reconstructed  phase  shift  (a)  and  iso-phase  contours  (b)  tank  pressure  of  10  psig 
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Figure  7-30  : Dark  central  ground  (a)  and  Schlieren  (b)  images  for  tank  pressure  of  10  psig 
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should  be  considered  only  as  first  order  estimates.  This  final  step  is  included  here  for  the 
sake  of  completeness  and  as  a demonstration  of  how  the  measured  data  would  be  used  in 
a fully  operational  tomographic  measurement  system.  The  convolutional  back-projection 
(CBP)  algorithm  was  used  due  to  its  speed.  The  reconstructed  phase  profiles  measured 
above  represent  a 97-by-97  sample  grid  with  a sample  spacing  of  0.47  mm.  In  terms  of 
the  tomographic  reconstruction,  this  represents  97  slices  of  the  test  region  with  each  slice 
oriented  normal  to  and  at  a different  position  on  the  presumed  axis  of  symmetry.  The 
CBP  algorithm  used  65  points,  centered  on  the  longitudinal  axis  of  the  flow,  and  the 
single  set  of  projection  data  was  duplicated  18  times  to  simulate  a 10  degree  projection 
angle  spacing.  Figure  7-3 1 shows  the  geometry  of  the  reconstructions  described  above. 

In  order  to  fully  reconstruct  the  density  profile  of  the  test  region,  97  separate 
tomographic  reconstructions  are  required  with  each  resulting  slice  representing  a 65 -by- 
65  point  cross  section  of  the  flow  region.  As  discussed  in  Chapter  3,  the  output  of  the 
tomographic  reconstruction  will  be  the  index  of  refraction  difference,  n-n^.  The  density 
ratio  can  be  calculated  using  the  Gladstone-Dale  relationship  as  follows. 

^.1-^  [7.1] 

Po 

Figures  7-32  And  7-33  are  plots  of  the  reconstructed  density,  centered  on  the  axis 
of  symmetry,  along  the  axis  of  flow  for  tank  pressures  of  30  and  20  psig.  These  were 
obtained  by  piecing  together  the  central  row  of  data  from  each  of  the  97  reconstructed 
slices.  The  tank  pressure  of  30  psig  represents  a supersonic  flow  condition  as  evidenced 
by  the  presence  of  shocks  at  the  exit  while  the  20  psig  tank  pressure  represents  a subsonic 
flow  condition.  Because  of  non-isentropic  flow  within  the  nozzle  it  is  not  possible  to 
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Figure  7-31  : Tomographic  reconstruction  geometry 


calculate  the  theoretical  density  ratio  at  the  nozzle  exit.  A more  theoretically  tractable 
test  setup  would  be  required  to  validate  the  reconstructed  density  ratio. 

Summary  and  Conclusions 

In  spite  of  the  constraints  imposed  by  the  choice  of  detector  array  and  COTS  lens 
arrays,  and  the  distortion  errors  which  resulted,  a two-dimensional  Shack-Hartmann 
wavefront  sensor  was  built  and  tested  with  some  success.  The  sensor  met  the  20  dB 
dynamic  range  specification  over  a 16  square  mm  portion  of  the  original  53  mm  test 
region.  The  results  of  the  free  air  jet  tests  demonstrate  the  measurement  capabilities  of 
the  system.  The  results  obtained  for  the  subsonic  flow  conditions  were  especially 
promising  since  this  is  an  area  where  other  flow  visualization  techniques  have  not 
performed  well.  Although  the  current  sensor  was  not  designed  for  this  measurement. 
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Figure  7-32  : Reconstructed  density  ratio  — (a)  and  iso-density  contours  (b)  tank  pressure  of  30  psig 
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Figure  7-33  : Reconstructed  density  ratio  — (a)  and  iso-density  contours  (b)  tank  pressure  of  20  psig 

\Po  / 


good  phase  data  was  obtained  utilizing  only  the  lower  portion  of  its  measurement 
dynamic  range. 


CHAPTER  8 

SUMMARY  AND  CONCLUSIONS 
General  Conclusions 

The  four  major  steps  required  to  carry  out  the  tomographic  reconstruction  of  a 
density  field  using  a wavefront  sensor  are  tomography,  phase  reconstruction,  spot 
deflection  measurement,  and  wavefront  sensor  design.  In  summarizing  this  work,  there 
are  a number  of  key  points  to  be  emphasized  from  each  of  these  four  subjects.  As  a 
whole  this  work  was  met  with  limited  success,  due  in  large  part  to  the  constraint  of  using 
an  off-the-shelf  lens  array.  In  spite  of  the  limitations,  useful  and  insightful  measurements 
were  made  of  a free  air  jet,  and  the  groundwork  has  been  laid  for  the  use  of  this  technique 
in  other  applications. 

Tomographic  Reconstruction 

The  key  point  to  be  made  concerning  the  tomographic  reconstruction  is  that  the 
results  are  dependent  upon  the  characteristics  of  the  measurement  being  made.  As  shown 
in  Chapter  3,  the  prediction  of  noise  variance  gives  a rough  estimate,  but  the  actual  noise 
performance  can  be  much  better  than  predicted  and  depends  upon  the  algorithm  used.  In 
this  work,  the  iterative  ART  algorithm  had  the  best  performance,  although  it  is  quite 
computationally  complex.  Also,  if  a high  resolution  reconstruction  is  desired  then  the 
number  of  projections  becomes  prohibitively  large  for  a high  speed  dynamic 
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measurement.  However,  it  has  been  shown  that  the  general  structural  details  can  be 
reconstructed  with  as  few  as  18  projections  spaced  at  10  degrees. 

Wavefront  Phase  Reconstruction 

A piecewise  least  squares  fit  utilizing  singular  value  decomposition  was  used  to 
reconstruct  phase  from  the  measured  phase  slopes.  As  discussed  in  Chapter  4,  this 
technique  had  marginal  performance  and  was  not  well  suited  to  applications  where  the 
measurement  contains  obscurations  or  sharp  discontinuities.  Further  work  is  needed  in 
this  area  to  develop  more  robust  algorithms. 

Spot  Deflection  Measurement 

At  the  heart  of  the  Shack-Hartmann  wavefront  sensor  is  the  process  of  accurately 
measuring  spot  position  and  deflection.  Errors  in  this  process  come  primarily  from 
detector  noise  and  from  the  algorithm  itself  With  any  centroiding  or  curve  fitting 
technique,  correction  of  the  algorithm  error  is  crucial.  The  4*’’  order  curve  fit  algorithm 
offers  the  best  performance  because  its  algorithm  error  is  zero  at  the  center  as  well  as  the 
edge  of  its  range.  In  choosing  a detector  array  and  lens  array,  the  camera  noise  must  be 
balanced  with  the  residual  algorithm  error  by  optimizing  the  spot  size  and  the  number  of 
detector  pixels  used  in  the  position  measurement  algorithm.  The  camera  noise  induced 
error  decreases  as  the  number  of  detector  pixels  is  decreased,  while  the  algorithm  error 
increases.  For  a given  number  of  utilized  detector  pixels,  the  camera  noise  induced  error 
has  a distinct  optimum  operating  point  in  terms  of  spot  size,  while  the  algorithm  error 
decreases  as  the  spot  size  is  increased.  For  the  Kodak  Megaplus  1 .4  camera  used  in  this 
work,  the  optimum  spot  size  was  around  4 detector  pixels  wide. 
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Wavefront  Sensor  Design 

It  is  difficult  to  outline  a generic  wavefront  sensor  design  procedure  because  of 
the  diverse  nature  of  the  measurement  applications  and  their  requirements.  Instead  the 
design  must  be  tailored  to  the  application  and  any  system  constraints.  The  relationships 
outlined  in  Chapter  6 provide  the  foundation  upon  which  any  design  can  be  carried  out. 
One  key  conclusion  which  should  be  emphasized  is  the  need  to  eliminate,  whenever 
possible,  any  optical  components  between  the  lens  array  and  detector  array.  As  seen  in 
Chapter  7 the  resulting  aberration  induced  errors  from  these  optics  can  be  very  large. 

Future  Research 

At  the  conclusion  of  this  work  there  are  several  areas  in  which  future  research  is 
warranted  and  should  lead  to  further  improvements  in  the  technique.  As  mentioned 
above,  the  phase  reconstruction  from  phase  slopes  is  in  need  of  further  improvement. 
Iterative  techniques  would  be  well  suited  to  this  application  and  should  be  explored  as 
one  possible  solution.  Beyond  this,  considerable  improvement  should  be  achieved  by 
using  a custom  lens  array  tailored  to  the  detector  array,  thus  eliminating  all  relay  optics 
between  the  two.  The  practical  aspects  of  mounting  and  configuring  the  lens  array 
detector  array  pairs  for  use  in  a full  tomographic  reconstruction  system  need  to  be 
examined. 

Characterization  of  the  system  was  accomplished  by  using  a linear  wavefront  tilt 
in  conjunction  with  an  independent  measurement  of  the  tilt  angle.  Although  this  was 
effective  for  this  proof  of  concept,  it  would  be  unwieldy  to  accomplish  in  a full 
tomographic  measurement  system.  A mechanical  tilting  of  the  sensor  itself,  as  discussed 
in  Chapter  7,  should  prove  to  be  a more  effective  means  of  calibrating  and  characterizing 
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the  wavefront  sensors.  In  addition,  a standard  calibration  test  object  with  a known  index 
profile  would  be  desirable  to  test  the  full  tomographic  reconstruction  system.  Further 
research  is  needed  to  examine  the  feasibility  of  developing  such  a standard. 

Beyond  these  specific  areas  of  improvement,  the  next  logical  step  is  to  develop  a 
full  tomographic  measurement  system  and  verify  that  valid  quantitative  data  can  be 
obtained.  The  quantitative  measurements  of  density  ratio  performed  in  this  work  are 
somewhat  tenuous  due  to  the  fact  that  only  a single  projection  was  available.  Before  an 
operational  system  could  be  feasible,  the  packaging  issues  discussed  above  need  to  be 
resolved.  Also,  the  cost  of  such  a system  could  be  cut  substantially  if  a laser  diode  could 
be  utilized  as  the  source.  The  power  budget  for  the  system  built  in  this  work  would  have 
allowed  the  use  of  a laser  diode  source. 

One  of  the  most  promising  results  to  come  from  this  work  was  the  ability  to 
obtain  useful  measurements  with  low  speed  flows  although  the  system  was  not  designed 
with  this  type  of  measurement  in  mind.  The  system  was  designed  to  provide  valid 
measurements  near  shock  boundaries  so  that  the  low  speed  flow  measurements  utilized 
only  the  lower  portion  of  the  systems’  full  dynamic  range.  This  is  an  application  for 
which  this  technique  is  well  suited.  By  utilizing  a camera  with  more  dynamic  range  than 
the  uncooled  8-bit  camera  used  in  this  work,  and  designing  the  sensor  for  higher 
sensitivity,  it  is  believed  that  this  type  of  system  could  prove  to  be  a very  useful  tool  in 
analyzing  low  speed  flows.  The  lack  of  shock  discontinuities  in  these  flow  measurements 
alleviates  many  of  the  difficulties  encountered  in  measuring  supersonic  flows. 
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APPENDIX 

MATLAB  SCRIPT  FILES 


SSOCFC.M  Hartmann  Array  Spot  Deflection  Calculations 

% ssqcfc.m 

% find  spot  centroid  shifts 

% Jay  Land  3/12/99 

clear 

close 


tic 

M=2; 

SPC=9; 

pw=SPC; 

lenap=188 ; 

pixwid=6 . 8 ; 

ref f ile= ' clnl28  ' ; 

testfile='cln058  ' ; 

% centroid  footprint  size=2M+l 
% lenslet  footprint  in  CCD  pixels 
% pad  all  boundaries  with  zeros  to  width  pw 
% lenslet  aperture  diameter  (microns) 

% CCD  pixel  width  (microns) 

disp( 'READING  REFERENCE  FILE') 

[imgdc,midpt] =readbmp3 (cat (2, ' c : \research\spots\dcdata\ ', ref file, ' . vl4 ' 
) ) ; 

disp( 'READING  TEST  FILE') 

imgtst=readbmp3b (midpt, cat (2, ' c: \research\spots\dcdata\ ' , testfile, ' . vl4 
' ) ) ; 


ilen=size (imgdc, 1) ; 
clen=size (imgdc, 2) ; 
N=floor (ilen/SPC) ; 
cptsx=zeros (N) ; 
cptsy=zeros (N) ; 
prowdc=zeros (N) ; 
pcoldc=zeros (N) ; 
prowtst=zeros (N) ; 
pcoltst=zeros (N) ; 
xcdc=zeros (N) ; 
ycdc=zeros (N) ; 
xctst=zeros (N) ; 
yctst=zeros (N) ; 
pvaldc=zeros (N) ; 
pvaltst=zeros (N) ; 
delxc  p=zeros(N); 
delyc  p=zeros(N); 
delxc  fp=zeros (N) ; 

% size  of  image  (rows) 

% size  of  image  (columns) 

% number  of  SPC  X SPC  sub-regions  sqrt 

delyc_fp=zeros (N) ; 

imgdcO=imgdc ( : , SPC+1 : SPC+ilen) ; 

imgtstO=imgtst ( : , SPC+1 : SPC+ilen) ; 

imgdc= [ zeros (pw, clen) ; imgdc; zeros (pw, clen) ] ; 

imgtst= [ zeros (pw, clen) ; imgtst ; zeros (pw, clen) ] 

row=pw+5;  % initial  dummy  value 
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column=pw-3 ; % Find  DC  image  peaks 

clc 

for  i=l:N  % row  index 

home 

dispCDC  PEAKS  % COMPLETE') 

100*i/N  % percent  complete 

if  i==l 
if  SPC==8 
row=13 ; 
col=5 ; 
end; 

if  SPC==9 
row=14 ; 
col=6; 
end; 
else 

row=prowdc (i-1, 1) +SPC; 
col=pcoldc (i-1, 1) -SPC; 
end; 

for  j=l:N  % column  index 

rorig=row- (floor (SPC/2) ) ; 

% origin=top  If  corner  of  new  search  region 
corig=col+ (floor (SPC/2) ) ; 

% column  origin  ( last  line  was  row  origin) 
subimg=imgdc (rorig : rorig+SPC-1 , corig : corig+SPC-1 ) ; 

% SPC  X SPC  search  sub-region 
[y,  ind] =max ( subimg,  [ ] , 2 ) ; 

% find  max  of  rows  & their  column  number  index 
[z,m] =max (y) ; % find  max  of  the  row  max's  and  row  index 

row=rorig+m-l ; % peak  value  row  index 

col=corig+ind (m) -1 ; % peak  value  column  index 

prowdc ( i , j ) =row; 
pcoldc ( i , j ) =col ; 
pvaldc ( i , j ) =subimg (m, ind (m) ) ; 
end; 

end 

%prowdc=prowdc-pw; 

%pcoldc=pcoldc-pw; 

% calculate  centroids  of  DC  image 
x=[-2;-l;0;l;2] ; 
for  i=l:N 
home 

disp('DC  CENTROIDS  % COMPLETE') 

100*i/N  % percent  complete 

for  j=l:N 

rs=prowdc (i, j ) -M; 
re=rs+ (2*M) ; 
cs=pcoldc ( i , j ) -M; 
ce=cs+ ( 2*M) ; 

csum=sum ( imgdc ( r s : re , cs  : ce ) , 1 ) ' ; 
r sum=sum ( imgdc ( r s : re , cs : ce ) , 2 ) ; 
a=polyf it (x, csum, 4 ) ; 
b=[4*a(l) ;3*a(2) ;2*a(3) ;a(4) ] ; 
rtsfit=roots (b) ; 

[pkabs, ind] =min (abs (rtsfit) ) ; 

xcdc (i, j ) =pcoldc (i, j ) +rtsfit (ind) -0.5; 

a=polyf it (x, rsum, 4 ) ; 
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b=[4*a(l) ;3*a(2)  ;2*a(3)  ;a(4)  ] ; 
rtsfit=roots (b) ; 

[pkabs, ind] =min (abs (rtsfit) ) ; 

ycdc (i, j ) =prowdc (i,j)+rtsfit(ind)-0.5; 

end; 

end; 

% Find  test  image  peak  locations 
for  i=l:N  % row  index 

home 

dispCTEST  PEAKS  % COMPLETE') 

100*i/N  % percent  complete 

for  j=l:N  % column  index 

rorig=prowdc ( i , j ) - ( floor { SPC/2 ) ) ; 

% origin=top  left  corner  of  new  search  region 
corig=pcoldc (i, j ) - (floor (SPC/2  ) ) ; 

% column  origin  ( last  line  was  row  origin) 
subimg=imgtst (rorig : rorig+SPC-1 , corig : corig+SPC-1 ); % SPC  X SPC 
search  sub-region 

[y,  ind] =max (subimg,  [] ,2) ; 

% find  max  of  rows  & their  column  number  index 
[ z, m] =max (y) ; % find  max  of  the  row  max's  and  row  index 

row=rorig+m-l ; % peak  value  row  index 

col=corig+ind (m) -1 ; % peak  value  column  index 

prowtst (i, j ) =row; 
pcoltst (i, j ) =col; 
pvaltst (i, j ) =subimg (m, ind (m) ) ; 
end; 

end 

%prowtst=prowtst-pw; 

%pcoltst=pcoltst-pw; 

% calculate  centroids  of  Test  image 
for  i=l:N 
home 

dispCTEST  CENTROIDS  % COMPLETE') 

100*i/N  % percent  complete 

for  j=l:N 

rs=prowtst (i, j ) -M; 
re=rs+ (2*M) ; 
cs=pcoltst (i, j ) -M; 
ce=cs+ (2*M) ; 

csum=sum (imgtst ( rs : re , cs : ce ) , 1 ) ' ; 
rsum=sum (imgtst (rs : re, cs : ce) ,2) ; 
a=polyf it (x, csum, 4 ) ; 
b=[4*a(l) ;3*a(2) ;2*a(3) ;a(4) ] ; 
rtsfit=roots (b) ; 

[pkabs, ind] =min (abs (rtsfit) ) ; 
xctst (i, j ) =pcoltst (i,j)+rtsfit(ind)-0.5; 
a=polyf it (x, rsum, 4 ) ; 
b=[4*a(l) ;3*a(2) ;2*a(3) ;a(4) ] ; 
rtsfit=roots (b) ; 

[pkabs, ind] =min (abs (rtsfit) ) ; 

yctst ( i, j ) =prowtst (i,j)+rtsfit (ind) -0.5; 

end; 

end; 

xcdcold=xcdc; 
ycdcold=ycdc; 
xctstold=xctst ; 
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yctstold=yctst ; 

% Correct  for  Centroiding  Algorithm  Error 

load  cfitcal;  % Load  calibration  Data 

(numoffs, centroid) 

realpos=linspace (0, 0 . 5, numof f s ) ' ; % Cal  data  actual  centroid  positions 


xcdcrel=xcdc-pcoldc+0 . 5 ; 
pixel  center 

ycdcrel=ycdc-prowdc+0 . 5 ; ' 


dc  field  x centroids  relative  to  peak 


dc  field  y centroids  rel  to  peak  pixel  center 
field  X cent  rel  to  pk  pix  center 
field  y cent  rel  to  pk  pix  center 


test 

test 


xctstrel=xctst-pcoltst+0 . 5 ; 
yctstrel=yctst-prowtst  + 0 . 5 ; 

CRSize=size (xcdcrel, 1) ; 

CCSize=size (xcdcrel, 2 ) ; 
for  i=l:CRSize 
home 

dispC  ERROR  CORRECTION  % COMPLETE') 

100*i/N 
for  j=l:CCSize 

xdcsrch=xcdcrel (i,  j ) ; 
ydcsrch=ycdcrel (i, j ) ; 
xtstsrch=xctstrel (i, j ) ; 
ytstsrch=yctstrel (i, j ) ; 

xcdcreln (i, j ) =sign (xdcsrch) *calcorr (sign (xdcsrch) *xdcsrch, numoffs, calcu 
rvf ) ; 


percent  complete 


ycdcreln (i, j ) =sign (ydcsrch) *calcorr (sign (ydcsrch) *ydcsrch, numoffs, calcu 
rvf)  ; 

xctstreln (i, j ) =sign (xtstsrch) *calcorr ( sign (xtstsrch) *xtstsrch, numoffs, c 
alcurvf ) ; 

yctstreln (i, j ) =sign (ytstsrch) *calcorr (sign (ytstsrch) *ytstsrch, numoffs, c 
alcurvf) ; 

end; 

end; 

xcdc=xcdcreln+pcoldc-0 . 5; 
ycdc=ycdcreln+prowdc-0 . 5; 
xctst=xctstreln+pcoltst-0 . 5; 
yctst=yctstreln+prowtst-0 . 5; 
delxc_p=xctst-xcdc; 
delyc_p=yctst-ycdc; 
sepldc=xcdc (1, N) -xcdc (1,1); 

sep2dc=xcdc (floor (N/2) , N) -xcdc (floor (N/2)  , 1)  ; 
sep3dc=xcdc (N, N) -xcdc (N, 1) ; 
sep4dc=ycdc (N, 1 ) -ycdc (1,1); 

sep5dc=ycdc (N, floor (N/2 ) ) -ycdc ( 1 , floor (N/2 ) ) ; 
sep6dc=ycdc (N, N) -ycdc ( 1 , N)  ; 

avgsepdc= ( sepldc+sep2dc+sep3dc+sep4dc+sep5dc+sep6dc) / 6; 
sepltst=xctst (1, N) -xctst (1,1); 

sep2tst=xctst ( floor (N/2 ) , N) -xctst ( floor (N/2 ) , 1 ) ; 
sep3tst=xctst (N, N) -xctst (N,  1 ) ; 
sep4tst=yctst (N, 1) -yctst (1,1)  ; 

sep5tst=yctst (N, floor (N/2) ) -yctst (1,  floor (N/2) ) ; 
sep6tst=yctst (N, N) -yctst { 1, N) ; 

avgseptst= ( sepltst+sep2tst+sep3tst+sep4tst+sep5tst+sep6tst ) /6; 

Mag_Cam= ( avgsepdc*pixwid) / ( (N-1) *lenap) ; 
delxc  fp=delxc_p*pixwid/Mag_Cam; 


X centroid  shift  in  CCD  pixels 
y centroid  shift  in  CCD  pixels 
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% X centroid  shift  in 
delyc_fp=delyc_p*pixwid/Mag_Cam; 

% y centroid  shift  in 
meanx=mean (mean (delxc_p) ) 
meany=mean (mean (delyc_p) ) 
stdx=std (delxc_p { : ) ) 
stdy=std (delyc_p ( : ) ) 
maxx=max (max (delxc_p) ) 
minx=min (min (delxc_p) ) 
maxy=max (max (delyc_p) ) 
miny=min (min (delyc_p) ) 
subplot (3,2, 1 ) 
mesh (delxc_p) 

axis([l  113  1 113  minx  maxx] ) 
title (sprintf (' %s  & %s  X-Shift 
subplot (3,2,2) 
mesh (delyc_p) 

axis([l  113  1 113  miny  maxy] ) 
title (sprintf (' %s  & %s  Y-Shift  (pixels 
subplot (3,2,3) 
contour (delxc_p) 
subplot (3,2,4) 
contour (delyc_p) 
subplot (3,2,5) 
hist (delxc_p ( : ) , 1000 ) 
grid  on 
subplot (3,2,  6) 
hist (delyc_p ( : 
grid  on 
toe/ 60 


at  lens  array  focal  plane 


at  lens  array  focal  plane 


(pixels) ',reffile,testfile) ) 


',reffile,testfile) ) 


) , 1000) 


function[img,midpt]=readbmp3 (vl4file)  ; 
fid=fopen (vl4file) ; 

[img, count] =fread ( f id,  7 6,  ' uintS  ' ) ; 

[img, count] =fread (fid, [1280 , inf] , ' uintS ' ) ; 
img=img ' ; 
f close ( fid) ; 

N=9; 

a=max (img (1 :N, : ) ) ; 
k=0; 

levelof f=0; 
amax=0 ; 
delamax=l ; 

while  delamax>0  & leveloff==0 
amax=max (a (1+ (k*N) : (k+1) *N) ) ; 
if  amax>50 

delamax=amax-amaxold; 
if  delamax<5  levelof f=l;  end; 
end; 

amaxold=amax ; 
k=k+l ; 
end; 

startpt= ( k-1 ) *N; 
a=fliplr (a) ; 
k=0; 

levelof f=0 ; 
amax=0 ; 
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delamax=l ; 

while  delamax>0  & leveloff==0 
amax=max (a ( 1+ ( k*N) : ( k+1 ) *N) ) ; 
if  amax>50 

delamax=ainax-amaxold; 
if  delamax<5  leveloff=l;  end; 
end; 

amaxold=amax; 

k=k+l; 

end; 

stoppt=size (a, 2) - ( (k-1) *N) ; 

midpt=round (startpt+ {0.5* (stoppt-startpt ) ) ) ; 
img=img ( : , midpt-512-N :midpt+511+N) ; 
return; 

function [img] =readbmp3b (midpt,  vl4file)  ; 
fid=f open (vl4 file ) ; 

[img, count ] =f read ( f id, 76, ' uintS ' ) ; 

[img, count] =f read (fid,  [1280, inf],  ' uint8 ' ) ; 
img=img ' ; 
f close ( fid)  ; 

N=9; 

img=img ( : , midpt-512-N : midpt +51 1+N) ; 
return; 

CERRCF2.M  4**^  Order  Curve  Fit  Algorithm  Error  Correction  Data  Generation 


clear 
close  all 

Pix  D=6.8;  % CCD  Pixel  Width  (urn) 

Sp_D=20;  % Measured  Spot  Null-Null  Width  (urn) 

Sp~P=Sp_D/Pix_D;  % Measured  Spot  Null-Null  Width  (CCD  Pixels) 

alpha=2/Sp_P; 
order=4 ; 

%alpha= . 5 ; 

x=linspace (-3, 3, 601) ' ; 
delx=x (2) -X (1) ; 
y=sincsq (x, alpha) ; 
plot (x, y) ; 
grid  on 

disp('Hit  Return  to  Continue') 
pause 

x=[-2,-l,0,l,2] ; 
for  os=-50:50 

offsetx(os+51)=os*delx; 

bin (1) =0.5* (2*sum(y(51-os:151-os) )-y(51-os)-y(151-os) ) ; 
bin (2 ) =0 . 5* (2*sum (y ( 151-os :251-os) )-y(151-os)-y(251-os) ) ; 
bin (3) =0 . 5* (2*sum (y (251-os : 351-os ) )-y(251-os)-y( 351-os ) ) ; 
bin ( 4 ) =0 . 5* ( 2*sum (y(351-os:451-os) )-y(351-os)-y(451-os) ) ; 
bin { 5 ) =0 . 5* ( 2*sum ( y ( 4 51-os : 551-os ) ) -y (4  51-os ) -y (551 -os ) ) ; 
a=polyf it (x, bin, order) ; 
for  i=l: order 

expn=order-i+l ; 
b (i) =expn*a  (i) ; 
end; 

rtsfit=roots (b) ; 
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[pkabs, ind] =min (abs (rtsfit) ) ; 
peak (os+51) =rtsfit (ind) ; 
end; 

err=peak-of fsetx; 
plot (of fsetx, err) 
grid  on 

offsetx2=offsetx(51:101) ; 
calcurvf=peak (51:101) ; 
numof f s=51 ; 

save ( ' of itcal ' , ' calcurvf ' , ' numof fs ' ) 

ZONAL  BUILD.M  SVD  Phase  Reconstruction  Algorithm 

function  phiM  = zonal_build  (delx,  dely,  L) 

% Mike  Gunia  11/99 
[m,n]  = size  (delx) ; 
o = m+ 1 ; 
of f set=0 ; 

sizel=floor (o/ (L-1) ) 
for  i=l:sizel-l 
for  j=l:sizel-l 
if  i==l  & j==l 

phi  = SVD_zonal  (delx ( (i-1) *L+1 : (i-1) *L+L-1, ( j-1) *L+1 : ( j-1) *L+L- 
1) , dely( (i-l)*L+l: (i-l)*L+L-l, (j-l)*L+l: (j-l)*L+L-l) ) ; 

(i-1) *L+1:  (i-1) *L+L 
(j-1) *L+1:  (j-1) *L+L 

phiM( (i-1) *L+1:  (i-1) *L+L,  (j-1) *L+1;  ( j-1) *L+L) =phi; 
end 

if  i==l  & j>l 

phi  = SVD_zonal  (delx((i-l)*L+l:(i-l)*L+L-l,(j-l)*L-(j-2):(j- 

1) *L-(j-2)+L-2),  dely ( (i-1) *L+1: (i-l)*L+L-l, ( j -1 ) *L- ( j -2 ) : (j-l)*L-(j- 

2) +L-2) ) ; 

(i-1) *L+1: (i-1) *L+L 
(j-l)*L-(j-2)  : (j-1) *L- (j-2)+L-l 

offset  = mean (phiM ( 1 : L, ( j -1 ) *L- ( j -2 ) ) -phi ( : , 1 ) ) ; 
phiM( (i-l)*L+l: (i-l)*L+L, ( j -1 ) *L- ( j -2 ) : ( j -1 ) *L- ( j -2 ) +L- 

1)  =phi+of fset; 

end 

if  j==l  & i>l 

phi  = SVD_zonal  (delx ( (i-1) *L- (i-2)  : (i-1) *L- (i-2) +L-2,  (j-1) *L- (j- 

2)  : (j-l)*L-(j-2)+L-2) , dely (( i-1 ) *L- ( i-2 ) : ( i-1 ) *L- ( i-2 ) +L-2 , (j-l)*L-(j- 
2)  : (j-l)*L-(j-2)+L-2) ) ; 

(i-1) *L- (i-2)  : (i-l)*L-(i-2)+L-l 
(j-1) *L+1: (j-1) *L+L 

offset  = mean (phiM (( i-1 ) *L- ( i-2 ),  j : j +L-1 ) -phi ( 1 ,:))  ; 
phiM( (i-1) *L- (i-2) : (i-1) *L- (i-2) +L-1, (j-l)*L+l: (j- 

1 )  *L+L) =phi+of f set ; 

end 

if  j>l  & i>l 

phi  = SVD_zonal  (delx ( (i-1 ) *L- (i-2 ):( i-1 ) *L- ( i-2 ) +L-2 , (j -1 ) *L- (j - 

2)  : ( j-1) *L- ( j-2) +L-2) , dely ( (i-1) *L- (i-2) : ( i-1 ) *L- ( i-2 ) +L-2 , ( j-1) *L- ( j- 
2) : (j-l)*L-(j-2)+L-2) ) ; 

(i-l)*L-i+2: (i-1) *L+L-l-i+2 
(j-l)*L-j+2: ( j-1) *L+L-l-j+2 

offset 1 = mean (phiM ((i-l)*L-i+2,  (j-l)*L-j+2:  (j-l)*L+L-l-j+2)- 
phi(l, :) ) ; 


off set 2 = mean (phiM {(i-l)*L-i+2: (i-l)*L+L-l-i+2, (j-l)*L-j+2) 
phi (:,!)) ; 

offset  = (of f set 1+of f set2 ) /2 ; 

phiM{ (i-1) *L-i+2: ( i-1 ) *L+L-l-i+2 , (j-1) *L-j+2: { j-1) *L+L-1- 
j +2 ) =phi+of f set ; 
end 

end 

end 

phiM=phiM(l: (i-l) *L+L-l-i+2, 1 : ( j -1 ) *L+L-1- j +2 ) ; 

[mnew, temp2] =size (phiM) ; 
if  (floor (mnew/2) ) ==mnew/2 
b = mnew/2-0.5; 
a = -b  ; 
else 

b = floor (mnew/2) ; 
a = -b ; 

end 

xl  = linspace  (a, b, mnew); 
yl  = xl; 

[X,  Y]  = meshgrid  (xl,yl); 

MAXX  = max(X( : ) ) ; 

MINX  = min(X( : ) ) ; 

MAXY  = max (Y ( : ) ) ; 

MINY  = min (Y ( : ) ) ; 

MINphiM  = min (min (phiM) ) 

MAXphiM  = max (max (phiM) ) 

phidiff=max  (max (phiM) ) -min  (min (phiM)); 
figure ( 1 ) 
mesh  (phiM) 

title  ('Zonal  Reconstruction  of  Flow  Field'); 
zlabel  ('Degrees'); 
xlabel ( ' (xl88um) ' ) ; 


function  phiM  = SVD_zonal  (xdev,  ydev) 

% Mike  Gunia  11/99 

%focal  = input  ('What  is  your  focal  length  in  microns?  : '); 
focal  = 10000; 

Sx  = xdev. /focal; 

Sy  = ydev. /focal; 
offset=0; 

[n,  tmpl]  = size (xdev);  %Size  of  actual  data  in  array 
m = n+1; 

% BUILD  A MATRIX  FOR  X SLOPES 

M = m*m;  %Length  to  include  all  elements  of  z data 

N = n*n;  %Length  to  include  all  elements  of  deviation  data 

Ax  = zeros (N, M) ; 

Ay  = zeros (N, M) ; 
for  i = 1:N 

if  mod ( i+of f set , m) ==0 
offset  = offset+1; 


end 

Ax (i,  i+offset) =-l;  % 
Ax (i, i+offset+1) =-l;  % 
Ax (i, i+offset+m) =1;  % 
Ax  ( i , i+of f set+m+1 ) =1 ; % 


% Ax  MATRIX  WILL  HAVE  THE  FORM 


-1 

-1 

0 

1 

1 

0 

0 

0 

0 

0 

-1 

-1 

0 

1 

1 

0 

0 

0 

0 

0 

0 

-1 

-1 

0 

1 

1 

0 

0 

0 

0 

0 

-1 

-1 

0 

1 

1 

Ay  MATRIX  WILL  HAVE  THE  FORM 
1-1  10  -1  100001 


Ay (i, i+offset)=-l; 
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Ay(i,i+offset+l)=l; 

% 1 

0 

-1 

1 

0 

-1 

1 

0 

0 

0 

Ay (i, i+offset+m) =-l; 

% 1 

0 

0 

0 

-1 

1 

0 

-1 

1 

0 

Ay (i, i+of f set+m+1 ) =1 ; 

% 1 

0 

0 

0 

0 

-1 

1 

0 

-1 

1 

end 

A = cat  ( 1 , Ax, Ay) ; 

A = A*0.5; 

%input  ( 'A  Matrix  complet; 

hit  any 

key 

to 

continue : 

' ) 

t 

%%%%%%%%% 

[u,d, v]=svd(A) ; % Do  SVD  to  get  inv[A]  %%%%%(75%  TIME  WEIGHT) %%%%% 

d; 

Dinv=zeros (M) ; % svd  inverse  for  d 

for  i=l:M 

if  d(i,i)<10"-6 

Dinv (i, i) =0; 

else 

Dinv (i, i ) =l/d { i, i ) ; 

end; 

end; 

u=u(:,l;M);  % u should  have  been  size (A)  ?? 

ut=u ' ; 

Ainv=v*Dinv*ut;  %%%%%  (25%  TIME  WEIGHT) %%%%% 

Sx  = (360/ (532*10^-7) )* (188*10^-4) * (tan(Sx) ) ; 

Sy  = (360/ (532*10''-7)  ) * (188*10-'-4)  * (tan(Sy)  ) ; 

%%%%%%%%%%%%% 

Sx  = Sx ( : ) ; 

Sy  = Sy ( : ) ; 

S = cat ( 1 , Sx, Sy) ; 

phi= (Ainv*S) ' ; % solve  for  xLS  - phase  values 

%timeSOL=toc 

%phi  = (A\S) ' ; 

phiM  = zeros  (m,m);  %phiM  is  the  MATRIX  loaded  version  of  wavefront 
data 

for  j = l:m 

phiM(j,:)  = phi ( j *m-m+l ; j *m)  ; 

end 

phiM  = phiM';  % ADD  OFFSET  TO  FINAL  PHI  VALUES 

CTFBP.M  Convolutional  Backproiection  Tomography  Algorithm 


% ctfbpb  Jay  Land  11/98 

% Filtered  Backpro j ection  Tomography  Algorithm 
% ESI4567  Project 


clear; 

tic; 

flp0=flops;  % 

%ftype=l;  % 

ftype=2;  % Shepp 

N=32; 

M=18; 

proj=zeros (4*N+1,M) ; 
f ilterh=zeros (4*N+1,  1) 
fildat=zeros (8*N+1,M) ; 
bakpro=zeros (8*N+1) ; 
rsize=4  *N+1 ; 
recon=zeros (rsize) ; 
rdat=zeros (2*N+1) ; 


Track  number  of  flops  required 

Ramachandran/Lakshminarayanan  direct  |f|  filter 
Logan  sinc(f/2fm)  windowed  filter  (2/pi  at  fmax) 
% 2N+1  points  used 

% Number  of  projection  angles  used 
% Initialize  projection  data  vector 
% Initialize  projection  filter  impulse  response 
% Initialize  filtered  projection  data  vector 
% Initialize  back-projected  data  vector 
% Size  of  reconstruction 
% Initialize  reconstruction  array 
% Truth  data  for  comparison 
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j=sqrt (-1) ; 
lam=632 . 8e-9; 
ko=2*pi/lam; 
no=l  ; 
nl=1.01; 
a=l/N; 
deln=nl-no; 
f lpl=f lops-f IpO ; 
filterh=ctf liter (N, ftype) ; 
flp2=flops-f IpO-f Ipl  ; 

[proj , rdat] =projg2 ( 4*N+1 , lam, M) 
flp3=f lops-f IpO-f Ipl- flp2  ; 
for  1=1 :M 

flldat ( : , 1) =a*conv (proj ( : , 1) , fllterh) ; 


wavelength  (meters) 

Free  space  wave  number 
base  Index 
peak  Index 

projection  data  sampling  Interval  (meters) 


% Filter  Impulse  response 
% gausslan  Index  profile 


Convolve  projection  data 


with  fllterh 
end; 

f lp4=f lops-f IpO-f Ipl- flp2-flp3; 
recon=bulldup (M, N, flldat ) ; 
flp5=flops-f IpO-f Ipl-f lp2-f Ip3-flp4  ; 

recon2= (pi/ (ko*M) ) *recon (N+1 : (3*N) +1, N+1 : (3*N) +1) ; % Scale  to  get  n-no 

reconstruction 

recon3= (pi/ (ko*M) ) *recon; 

mesh (recon2) ; 

dn=recon2 (N+1, N+1)  % Find  radial  peak  value 

err=rdat-recon2;  % Errors  In  reconstruction 

peak_error=max (max (err) ) % Determine  peak  error 

rmserr=sqrt  ( ( sum  ( sum  ( err  . ''2 ) ) ) / ( 2 *N+1 ) ^2 ) % RMS  Error 

rec2= (recon2-mln (min (recon2 ) ) ) *255 /max (max (recon2) ) ; 
rec3= (recon3-mln (min (recon3) ) ) *255/max (max (recon3) ) ; 
dlsp (sprlntf (' Elapsed  Time  (min)  %1 . 2f ' , toc/60) ) ; 
f lp6=f  lops-f  IpO-f  Ipl-f  lp2-f  lp3-f  lp4-f  lp5  ; 
flptot=f lops-f IpO ; 


function  fllterh=ctfllter (N, ftype) ; 
f llterh=zeros (4*N+1,  1)  ; 

a=l/N;  % Sample  spacing 

for  1=1:4*N+1;  % Define  filter 

n=l-l-2*N; 

If  ftype==l  % Ramachandran/Lakshmlnarayanan  Filter 

If  n==0  f llterh (1 ) =1/ ( 4*a^2 ) ; else 

If  (n/2)==round(n/2)  f llterh ( 1 ) =0 ; else  fllterh(l)=- 
1/ ( (n"2)  * (pl^2)  * (a''2)  ) ; end; 
end; 
end; 

If  ftype==2  % Shepp/Logan  Filter 

fllterh (1)= (-2/ ( (pl"2) * (a^2) ) ) / ( (4*n^2)-l) ; 
end; 
end; 

function  [proj , rdat] =projg2 (D, lam, M) ; 

Ql=2;  % Gaussian  width  Roll  Off  Factor  In  X direction 

q2=4;  % Gaussian  width  Roll  Off  Factor  In  Y direction 

alpha=2.5;  % Frequency  scale  factor 

N=(D-1) /4; 
proj=zeros (D,M) ; 
rdat=zeros (2*N+1)  ; 


% Initialize  projection  data  vector 
% Truth  data  for  comparison 
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ko=2*pi/lam; 
no=l ; 
nl=1.01; 
deln=nl-no; 
for  p=l:M 

g ainina =p  * p i / M ; 
cganuna=cos  (gamma)  ; 
sgamma=sin (gamma) ; 
for  i=l:4*N+l; 
n=i-l-2*N; 

for  k=l:2*N+l; 


% Free  space  wave  number 
% base  index 
% peak  index 


% Projection  Angle 


m=k-l-N; 

x=m/N; 

y=n/N; 

xr=x*cgamma+y*sganraia;  % Rotated  coordinates 

yr=y* cgamma-x*s gamma; 
index (k) =no+deln*exp (- 

0.5* ( (Ql*xr ) ^2+ (Q2*yr ) ^2 ) ) *cos (alpha*pi*xr ) ^2; 
end; 

proj (i, p) =CatQNC (index, -1, 1, 3, N+1) ; % Integration 

end; 
end; 

pro j =pro j * ko ; % Scale  for  phase 

for  m=l:2*N+l  % Compute  Truth  data 

y= (-N+m-1) /N; 
for  n=l:2*N+l 


x= (-N+n-1) /N; 

rdat  (m,n)=deln*exp(-0.5*  ( (Ql*x)  ^2+  (Q2*y)  "2)  ) *cos  (alpha*pi*x)  ''2; 
end; 
end; 
return; 


function  recon=builup (M, N, f ildat ) ; 

recon=zeros (4*N+1) ; % Initialize  reconstruction  array 

[x, y] =meshgrid (-4*N: 1 : 4*N) ; % Define  unrotated  data  grid 

[xr, yr] =meshgrid (-2*N: 1 : 2*N) ; % Define  unrotated  reconstruction  grid 

for  p=l:M  % Backproj ection  Summation 

gamma=p*pi/M;  % Projection  rotation  angle 

xm=xr*cos (gamma) -yr*sin (gamma) ; % Rotate  reconstruction  grid 

ym=xr*sin (gamma) +yr*cos (gamma) ; 

for  i=l:8*N+l  % Backproject  the  filtered  data 

bakpro ( : , i) =fildat ( : , p) ; 
end; 

f lpinit=f lops  ; 

recon=recon+interp2 (x, y, bakpro, xm, ym, ' cubic ' ) ; 

% Interpolate  onto  recon  grid 

flp (p) =flops-flpinit; 

mesh (recon (33 : 98, 33 : 98 ) /p)  % Plot  intermediate  result 

pause (.01)  % Delay  for  plotting 

end; 

pause ( 2 ) 

function  numi  = CatQNC ( fname, a, b, m, n) 

% 

% Pre: 

% fname  string  that  names  an  available  function  of  the 


% 

% 

% 


form  f(x)  that  is  defined  on  [a,b].  f should 
return  a column  vector  if  x is  a column  vector. 
a,b  real  scalars 
% m integer  that  satisfies  2 <=  m <=11 

% n positive  integer 


% 

% 

% 

% 

% 


Post : 

numi  the  composite  m-point  Newton-Cotes  approximation  of  the 
integral  of  f(x)  from  a to  b.  The  rule  is  applied  on 
each  of  n equal  subintervals  of  [a,b]. 

Delta  = (b-a) / (n-1) ; 
h = Delta/ (m-1) ; 

%x  = a+h* (0 : (n* (m-1) ) ) ' ; 
w = WNC (m) ; 

X = linspace (a, b, 2*n-l ) ' ; 
f=fname ' -1; 

% f = feval ( fname, x) ; 
numI  = 0 ; 
first  = 1; 
last  = m; 
for  i=l:n-l 

%Add  in  the  inner  product  for  the  i-th  subintegral, 
numi  = numi  + w' *f (first : last) ; 
first  = last; 
last  = last+m-1; 

end 

numi  = Delta*numl; 


function  w=WNC (m) ; 

% 

% Pre : 

% m integer  that  satisfies  2 <=  m <=  II 
% 

% Post : 

% w column  m-vector  consisting  of  the  weights  for  the 
% the  m-point  Newton-Cotes  rule. 

% 

if  m==2 

w=[l  1] ' /2; 
elseif  m==3 

w=[l  4 1] ' /6; 
elseif  m==4 

w=[l  3 3 1] ' /8; 
elseif  m==5 

w=[7  32  12  32  7] ' /90; 
elseif  m==6 

w=[19  75  50  50  75  19]  ' /288; 
elseif  m==7 

w=[41  216  27  272  27  216  41]'/840; 
elseif  m==8 

w=[751  3577  1323  2989  2989  1323  3577  751]'/17280; 
elseif  m==9 

w=[989  5888  -928  10496  -4540  10496  -928  5888  989]'/28350; 
elseif  m==10 

w=[2857  15741  1080  19344  5778  5778  19344  1080  15741  2857]'/89600 
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else 

w=[16067  106300  -48525  272400  -260550  427368  -260550  272400  -48525 
106300  16067] ' Z598752; 
end; 


Iterative  ART  Tomography  Algorithm 


% ctartl2  Jay  Land  5/99 
% Iterative  ART  Algorithm 
% Full  ART  Iterative  Reconstruction 
% Program  Calculates  the  weighting  coefficients 
% based  on  the  area  of  the  intersection  between 
% each  reconstruction  pixel  and  the  path  of  integration 
% which  is  one  pixel  wide.  The  weight  for  reconstruction 
% pixel  i,j  is  the  ratio  of  the  area  intersected  by  that 
% pixel  to  the  total  area  of  the  integration  path. 


clear; 
load  Pdata 
load  projfile 
tic; 

N=56; 

M=18; 

proj=zeros (4*N+1,M) ; 
rsize=2*N+l ; 
recon=zeros (rsize) ; 
recon2=zeros (rsize) ; 
rdat=zeros (rsize) ; 
errhist=zeros (30, 1) ; 
iteration 

rechis=zeros (rsize, rsize, 30) ; 

lam=532e-9; 

ko=2*pi/lam; 

no=l ; 

nl=1.01; 

a=l/N; 

interval  (meters) 

stdproj=0; 

for  i=l:M 


% 2N+1  points  used 

% Number  of  projection  angles  used 
% Initialize  projection  data  vector 
% Size  of  reconstruction 
% Initialize  reconstruction  array 

% Truth  data  for  comparison 
% history  of  rms  error  after  each 


% wavelength  (meters) 

% Free  space  wave  number 
% base  index 
% peak  index 

% projection  data  sampling 


proj ( : , i)  = [zeros (N, 1) ; Proj_Dat ; zeros (N,  1) ] +stdpro j *randn (4*N+1, 1) ; 
end; 

Lin_Ind= [Lin_Ind ( 1 ) ;Lin_Ind]  ; 
for  i=l:2*N+l 


ilin=i-N-l ; 
for  j=l:2*N+l 
j lin=j -N-1 ; 

indval=sqrt  ( ilin''2  + j lin"'2 ) ; 
if  indvaKN 

indlow=floor (indval)  ; 
indhigh=ceil (indval)  ; 

indslope=Lin_Ind ( 1+indhigh) -Lin_Ind ( 1+indlow) ; 
rdat (i, j ) =Lin_Ind ( 1+indlow) +indslope* (indval-indlow) ; 
else 

rdat ( i , j ) =Lin_Ind ( 1+N ) ; 
end; 
end; 
end; 
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rdat=Ind_o-rdat  ; 
itnum=l ; 
rmserr=l; 
rmse_old=2; 
keep_on=l; 
while  keep_on==l 
for  m=l : M 

gamma=in*pi/M  % Projection  Angle 

itnum 

if  gamma==0  pcase=0;  end; 

if  0<gamina  & gamma<=pi/4  pcase=l;  end; 

if  pi/4<gamma  & gairtma<pi/2  pcase=2;  end; 

if  gamma==pi/2  pcase=3;  end; 

if  pi/2<gamma  & gamma<=3*pi/4  pcase=4;  end; 

if  3*pi/4<ganuna  & gamma<pi  pcase=5;  end; 

if  ganuna==pi  pcase=6;  end; 

switch  pease 

case  0 

theta=0 ; 
for  n=l:4*N+l 

w=zeros (2*N+1) ; 
for  i=l:2*N+l 
for  j=l:2*N+l 

yc=abs { j-n+N) ; 
if  yc<2 

w (i, j ) =Proj_Area (theta,  yc)  ; 
end; 
end; 
end; 

sumw=sum (sum (w) ) ; 
if  sumw>0 

Pcalc=sum ( sum ( w . * recon) ) ; 

Perror=proj (n,m)-Pcalc; 
recon=recon+Perror*w/ sumw; 
end; 
end; 
case  1 

theta=gamma; 
for  n=l : 4 *N+1 

w=zeros (2*N+1) ; 
for  i=l:2*N+l 


for  j=l:2*N+l 

yc=abs ( ( ( ( (n-2*N-l) /sin (theta) ) -i+N+1) * tan (theta) ) - 

j+N+1) ; 

if  yc<2 

w (i, j ) =Proj_Area (theta,  yc) ; 
end; 
end; 
end; 

sumw=sum ( sum ( w) ) ; 
if  sumw>0 

Pcalc=sum ( sum ( w . * recon) ) ; 

Perror=proj (n,m)-Pcalc; 
recon=recon+Perror*w/ sumw; 
end; 
end; 
case  2 


theta= (pi/2 ) -gamma; 
for  n=l:4*N+l 

w=zeros (2*N+1) ; 
for  i=l:2*N+l 
for  j=l:2*N+l 

yc=abs ( ( (N+l-j ) /tan (gamma) ) + ( (n-2*N-l) / sin (gamma) ) 

+N+1) ; 

if  yc<2 

w (i, j ) =Proj_Area (theta, yc) ; 
end; 
end; 
end; 

sumw=sum ( sum ( w) ) ; 
if  sumw>0 

Pcalc=sum(sum(w. *recon) ) ; 

Perror=proj (n,m)-Pcalc; 
recon=recon+Perror*w/ sumw; 
end; 
end; 
case  3 

theta=0 ; 
for  n=l:4*N+l 

w=zeros (2*N+1) ; 
for  i=l:2*N+l 
for  j=l:2*N+l 

yc=abs (i-ntN) ; 
if  yc<2 

w (i,  j ) =Proj_Area (theta, yc) ; 
end; 
end; 
end; 

sumw=sum (sum (w) ) ; 
if  sumw>0 

Pcalc=sum (sum (w. *recon) ) ; 

Perror=proj (n,m) -Pcalc; 
recon=recon+Perror*w/sumw; 
end; 
end; 
case  4 

theta=gamma- (pi/2) ; 
for  n=l:4*N+l 

w=zeros (2*N+1) ; 
for  i=l:2*N+l 
for  j=l:2*N+l 

yc=abs ( ( ( j-N-1) *tan (theta) ) + ( (n-2*N-l) /cos (theta) ) 

+N+1) ; 

if  yc<2 

w (i, j ) =Proj_Area (theta,  yc)  ; 
end; 
end; 
end; 

sumw=sum ( sum ( w) ) ; 
if  sumw>0 

Pcalc=sum ( sum ( w . * recon) ) ; 

Perror=proj (n,m) -Pcalc; 
recon=recon+Perror*w/sumw; 
end; 
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end; 
case  5 

t he t a=p i - gamma ; 
for  n=l:4*N+l 

w=zeros (2*N+1) ; 
for  i=l:2*N+l 
for  j=l:2*N+l 

yc=abs ( ( (i-N-1) /tan ( (pi/2) -theta) ) - { (n-2*N- 
1) / (tan ( (pi/2) -theta) *cos ( (pi/2) -theta) ) ) -j+N+1) ; 
if  yc<2 

w (i, j ) =Proj_Area (theta,  yc)  ; 
end; 
end; 
end; 

sumw=sum ( sum (w) ) ; 
if  sumw>0 

Pcalc=sum ( sum (w. *recon) ) ; 

Perror=proj (n,m)-Pcalc; 
recon=recon+Perror*w/sumw; 
end; 
end; 
case  6 

theta=0 ; 
for  n=l:4*N+l 

w=zeros (2*N+1) ; 
for  i=l:2*N+l 
for  j=l:2*N+l 

yc=abs (3*N+2-n-j ) ; 
if  yc<2 

w (i, j ) =Proj_Area (theta,  yc)  ; 
end; 
end; 
end; 

sumw=sum ( sum (w) ) ; 
if  sumw>0 

Pcalc=sum ( sum ( w . *recon) ) ; 

Perror=proj (n,m)-Pcalc; 
recon=recon+Perror*w/ sumw; 
end; 
end; 
end; 
end; 

if  itnum>l  recold=recon2 ; end; 
recon2= ( 1/ko/a) * recon; 
rechis  itnum)  =recon2 ; 

mesh (recon2 ) ; 
pause ( 2 ) 

dn=recon2 (N+1, N+1)  % Find  radial  peak  value 

err=abs (rdat ' -recon2) ; % Errors  in  reconstruction 

err=err { 3 : 2*N-1 , 3 : 2*N-1 ) ; 

peak_error=max (max (err ) ) % Determine  peak  error 

rmse_old=rmserr; 

rmserr=sqrt  ( ( sum  ( sum  (err  . ''2 ) ) ) / ( 2*N+1 ) "'2 ) % RMS  Error 

errhist (itnum) =rmserr; 

rec2= (recon2-min (min (recon2) ) ) *255/max (max (recon2) ) ; 
disp (sprintf (' Elapsed  Time  (min)  % 1 . 2f ' , toe/ 60 ) ) ; 
if  itnum>3 
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if  rmserr<rmse_old 
keep_on=l ; 
else 

keep_on=0; 

end; 

else 

keep_on=l; 

end; 

itnum=itnum+l ; 
pause ( 2 ) 

end 

recon2=recold; 
mesh (recon2) ; 

dn=recon2 (N+1, N+1)  % Find  radial  peak  value 

err=abs (rdat ' -recon2 ) ; 

% Errors  in  reconstruction 
err=err ( 3 : 2*N-1 , 3 : 2*N-1 ) ; 

peak_error=max (max (err) ) % Determine  peak 

rmserr=sqrt  { (sum(sum{err . ''2)  ) ) / (2*N+1)  ^'2)  % RMS  Error 

itnum-1 


error 
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